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Abstract 


The  sum-difference  surface  current  formulation  is  introduced  to  treat  electromagnetic 
boundary  value  problems  when  anisotropic  impedances  are  specified  on  both  sides  of  a 
surface.  It  can  also  be  applied  to  impedance  coated  bodies.  This  formulation  preserves  the 
duality  nature  of  Maxwell  equations  and  carries  it  over  into  the  algebraic  form  of  the 
integrodifferential  operators  in  the  equations  for  surface  currents.  Since  a  90°  rotation  is 
equivalent  to  undergoing  a  duality  transform  for  an  incident  plane  wave,  this  particular 
symmetry  in  the  algebraic  form  of  the  operators  leads  to  sufficient  conditions  under  which 
the  on-axis  backscattering  of  an  anisotropic  impedance  coated  scatterer  having  a  90° 
rotational  symmetry  is  eliminated.  The  sum-difference  formulation  is  utilized  for  solving  the 
problem  of  electromagnetic  scattering  from  an  anisotropically  impedance  coated  tubular 
cylinder  of  finite  length.  The  solution  has  been  coded  in  FORTRAN  and  tested.  Some 
interesting  results  are  presented  and  discussed. 
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I.  INTRODUCTION 

Sometime  ago  the  question  was  raised:  "For  electromagnetic  boundary  value 
problems  with  specified  surface  impedances,  how  can  one  go  from  a  non-perfectly 
conducting  surface  on  which  both  the  electric  and  the  magnetic  equivalent  surface  currents 
are  to  be  found,  to  a  perfectly  conducting  surface  on  which  the  number  of  unknowns  is 
halved  [1]?"  The  answer  to  this  question  turns  out  to  be  one  of  algebra.  It  is  well  known  that 
the  impedance  specified  on  the  surface  of  a  body  separates  its  interior  completely  from  its 
exterior.  Therefore  an  impedance  coated  body  can  always  be  considered  as  a  hollow  volume 
enclosed  by  an  infinitesimally  thin  shell  with  surface  impedances  specified  both  on  the  inside 
and  the  outside  of  the  shell.  The  inside  and  the  outside  of  the  body  can  be  considered  as 
constituted  of  the  same  medium  and  the  impressed  electromagnetic  excitation  can  be  treated 
as  continuous  across  the  shell.  On  the  outside  surface,  there  are  the  equivalent  total  electric 

current  K +  and  total  magnetic  current  L + ;  on  the  inside  surface,  there  are  K  and  L  .  For 
an  exterior  problem,  only  K +  and  L +  need  to  be  found;  for  an  interior  problem,  only  K    and 

L    are  necessary.  A  single  formulation  for  solving  both  types  of  problems  would  appear  to 

require  finding  all  inside  and  outside  currents  therefore  doubling  the  amount  of  work,  but  it 
turns  out  not  to  be  the  case  because  some  of  the  currents  are  linear  combinations  of  others. 
Furthermore,  this  formulation  holds  the  key  to  answering  the  question  posed  above. 

Since  the  shell  is  infinitesimally  thin,  from  Maxwell  equations  the  radiation  to  the 
outside  and  to  the  inside  of  the  shell  can  both  be  given  in  terms  of  integrodifferential 

operators  on  the  sum  currents  K  =  K    +  K    and  L  =  L+  L   .  Note  that  the  outside  currents 

will  not  contribute  to  the  radiation  in  the  interior  while  the  inside  currents  will  not  contribute 
to  the  radiation  to  the  exterior  of  the  shell.  For  simplicity  in  the  description,  we  consider  the 
exterior  problem  of  electromagnetic  scattering.  By  definition,  the  radiation  is  the  difference 
between  the  total  field  and  the  incident  field.  On  the  surfaces  of  the  shell,  this  definition  links 

the  incident  E  and  H  fields  and  the  difference  currents  K    -  K    and  L-L    to  the  sum 

currents.lt  is  therefore  natural  to  treat  the  difference  currents  and  the  sum  currents  as  the  four 


1 


unknowns  to  be  solved  instead  of  the  inside  and  the  outside  currents. 

The  surface  impedance  on  the  outside  surface  of  the  shell  normalized  to  the  medium 
is  denoted  by  Z+  and  that  on  the  inside  surface  is  Z\  They  can  be  tensors  if  the  impedances 
are  anisotropic  and  may  vary  from  point  to  point.  By  forming  the  sum  impedance 
Z  =  (Z*  +Z")/2  and  the  difference  impedance  A  =  (Z"  -Z")/2,the  impedance  boundary 
conditions  provide  a  set  of  relations  between  the  difference  and  the  sum  currents.  It  turns  out 
that  the  rank  of  Z  determines  how  the  unknown  surface  currents  are  solved.  If  Z  is  invertible, 
then  the  difference  currents  are  linear  combinations  of  the  sum  currents  so  that  only  the 
integrodifferential  equations  of  the  sum  currents  have  to  be  solved.  There  are  only  two 
unknowns  to  be  solved  for  both  the  exterior  and  the  interior  problems  under  this  sum- 
difference  current  formulation.  If  Z  is  rank  0,  then  Z  =  0;  the  impedance  boundary  condition 

requires  that  L  be  proportional  to  a  90°  rotation  of  A^.  The  difference  electric  current  is 

obtained  from  K  after  the  integrodifferential  equation  on  K  is  solved.  This  situation  includes 

the  case  where  Z +  =  Z "  =0  when  the  surface  is  perfectly  conducting,  thus  the  result  answers 
the  question  about  the  transition  of  the  equations  from  a  problem  of  two  variables  to  one 
which  has  only  a  single  variable. 

Instead  of  dealing  with  an  impedance  coated  body,  this  thesis  presents  the  sum- 
difference  currents  formulation  of  electromagnetic  boundary  value  problems  for  the 
scattering  of  an  infinitesimally  thin  surface  for  which  both  the  inside  and  the  outside  currents 
are  true  unknowns  to  be  found.  Extension  of  this  formulation  to  impedance  coated  bodies  is 
then  discussed. 

This  formulation  preserves  the  duality  nature  of  Maxwell  equations  and  carries  it  over 
into  an  explicit  specific  algebraic  form  of  the  integrodifferential  operators  in  the  equations 
for  the  sum  currents.  Since,  for  a  plane  wave,  a  90°  rotation  is  equivalent  to  undergoing  a 
duality  transform,  this  explicit  symmetry  in  the  algebraic  form  of  the  operators  enables  us  to 
deduce  sufficient  conditions  for  eliminating  the  on-axis  backscattering  of  an  anisotropic 
impedance  coated  scatterer  having  a  90°  rotational  symmetry. 

The  sum-difference  currents  formulation  is  utilized  for  solving  the  problem  of 


electromagnetic  scattering  from  an  anisotropically  impedance  coated  tubular  cylinder  of 
finite  length.  The  solution  has  been  coded  in  FORTRAN  and  tested.  Some  interesting  results 
are  presented  and  discussed. 

In  this  thesis,  the  time  dependence  e "'"'  is  used.  E  represents  the  electric  field 

intensity  divided  by  the  intrinsic  impedance  of  the  medium,  \J\ile ;  therefore  it  has  the  same 

unit  as  H  in  amperes  per  meter.  So  are  the  electric  and  magnetic  equivalent  surface  currents  K 

and  L. 


II.  THE  SUM-DIFFERENCE  SURFACE  CURRENT  FORMULATION  OF 
ELECTROMAGNETIC  BOUNDARY- VALUE  PROBLEMS 

A.         STRATTON-CHU  FIELD  FORMULATION  AND  RADIATION 

On  an  orientable,  piecewise  regular  surface,  whether  open  or  closed,  having  a  surface 
electric  current  density  K  and  a  surface  magnetic  current  density  L,  we  define  the  Stratton- 
Chu  E-field  formula  as: 

EsJr,S,K,L)    =  £!   f^)G(f-r>rfv  f  K(ro)  -V  G(r-ro)dao 

-  ^V  x  [L(ro)(Kr-ro)dao 


. e<*|f"fJ 

where  f  is  a  point  which  is  not  on  S,  k  =  (dj\io£o  and  G(r-ro)  =  —  ,  then  the 


k\r-r 


Stratton-Chu  H-field  formula  can  be  defined  as: 
HsJF,S,K,L)  =  E(r,S,L,-K) 


Ol 


AV  x  fi(VGp.rjdam  *  g  /s^)G(?-?A       (2.2) 

"^v//(r°)'v''G(?"F">'ia° 


Note  that  if  5  is  a  closed  surface  and  AT  and  L  are  the  actual  total  equivalent  surface  currents 
on  S,  then  Es_c  and  Hs_c  are  respectively  the  E  and  H  fields  at  r  due  to  all  sources  inside 
S  if  r  is  located  outside  S  and  vice  versa.  This  is  a  direct  consequence  of  Maxwell's 
equations  [2]  and  under  this  circumstance,  the  Stratton-Chu  formulae  are  equivalent  to 
Maxwell's  equations.  On  the  other  hand,  unlike  the  Poynting  theorem,  the  Stratton-Chu  field 
formulae  over  an  open  surface  S  are  but  integrodifferential  operators  on  the  tangential  vector 


fields  K  and  L  over  S  without  any  special  physical  meaning  attached. 

To  introduce  equivalent  surface  currents  on  S,  the  direction  of  the  unit  normal  vector  n 
on  the  surface  has  to  be  chosen.  Adopting  the  terminology  for  a  closed  surface,  we  can  assign 
one  side  of  any  orientable  surface  S  as  the  "outside  surface"  S+,  albeit  somewhat  arbitrarily 
if  S  is  not  closed.  The  outward  normal  n  +  is  the  unit  normal  pointing  out  of  this  side  of  S 

and  for  simplicity,  we  call  n  +  the  normal  of  S  and  denote  it  by  n.  The  other  side  of  the 
surface  S  is  the  "inside  surface"  S~.  At  any  point  r  on  S,  the  inward  normal  n  ~  is  -n.  As  a 

convention,  the  fields  and  surface  currents  on  S  +  and  S '  always  carry  the  corresponding 
superscripts  (Figure  2-1). 


Figure  2-1.  Outside  and  Inside  Surfaces,  Normals  and  the  Equivalent  Currents. 


Each  of  the  total  surface  currents  K+  and  L*  on  5*  consists  of  two  parts:  the  incident 


current  (with  the  additional  superscript  "inc")  and  the  scattering  current  (with  the  additional 
superscript  "sc")  corresponding  to  the  incident  field  and  the  scattered  field  on  the  particular 
side  of  the  surface  S: 


K     =  n~  x  H     =  n~  x  H       +  n~  x  H 


±,sc 


(2-3) 


=  K*-™  +    K*sc 


L     =  £    x  n~  =  E       x  n~  +  E        x  n 


=  L         +    L 


(2-4) 


Note  that  S  is  infinitesimally  thin,  hence  H +'  c  =  H  '  lc  =  H  mc  and  E +'  Ic  =  E  '  lc  =  E  mc  on 

S  so  that  jf  +'OTC  =  -  K  'inc  and  L  +'inc  =  -  L  'inc .  Therefore  the  sum  currents  K  and  L  on  S 
defined  below  are  also  the  corresponding  sums  of  the  scattering  currents  only: 


K  =  K+  +  K~  =  K+'sc  +  K~'sc 


f*  f+  f*~  f*+>.sc  f-,sc 

L  =  L     +  L      =  L        +  L 


(2-5) 


Since  the  Stratton-Chu  field  formulae  are  linear  operators  on  the  surface  currents,  the 
radiated  fields  from  surface  currents  on  S  are  determined  by  the  sum  currents  only: 

Es\r)  =  EsJr,S\K\L+)  +  Es_c(r,S-,K'X')  =  EsJr,S,K,L)         _ 

(2-6) 
#  "(r)  -  Hs_c(r,S,K,L)  =  Es_c(r,S,L,-K) 


B.         CONDITION    ON    THE    CURRENTS    IMPOSED    BY    MAXWELL'S 
EQUATIONS 

As  r  approaches  r*  on  5*,  the  tangential  components  (denoted  by  the  subscript 


"tan")  of  Eq.  (2-6)  provide  four  equations  relating  the  incident  fields  and  the  total  currents 
on  both  sides  of  S  through  the  fact  that  the  incident  field  is  the  difference  between  the  total 
and  the  scattered  field: 

C  =  ^  -  Es_cm(r\SXL)  (2-7) 


C  =  Cn  ~  ES.C^-,S,K,L)  (2-8) 


HZ  =H+    -H(r,S,K,L)  (2-9) 


CC  =C  -  Hs_c,Jir-,S,K,L)  (2-10) 


These  four  equations  are  not  independent  of  each  other.  Because 
n*  x  (E(r*)  x  **)  =  gj?*)  =  n*  x  L± 
ft*  x  (H(r*)  x  n*)  =  Hm(r*)  =  -n*  x  K* 


(2-11) 


the  difference  between  Eqs.  (2-7)  and  (2-8)  trivially  confirms  the  definition  of  the  sum 
equivalent  magnetic  current  while  the  difference  between  Eqs.  (2-9)  and  (2-10)  confirms  the 
definition  of  the  sum  equivalent  electric  current  as  both  can  be  deduced  directly  from 
Maxwell's  equations.  We  choose  to  use  the  sum  of  Eqs.  (2-7)  and  (2-8)  and  that  of  Eqs.  (2-9) 
and  (2-10)  as  the  two  independent  linear  combinations  out  of  Eqs.  (2-7)  through  (2-10)  to 
link  the  incident  fields  to  the  total  surface  currents  on  5*  as  dictated  by  Maxwell's  equations: 


IE, 


tan 


2H, 


tan 


=  A  x  (L+  -  L~)  -  {EsJr\S,K,L)  +  E,_c(r,S,£L)} 

=  w  x  (L+  -  L")  +  MK  -  NL 

=  -  nx(K+  -  K')  -  {HsJr\S,K,L)  +  tf  _c.(r,S,£L)} 

=  -  n  x  (K*  -  K~)  +  NK  +  ML 


(2-12) 


where  M  and  ,/V  are  linear  integrodifferential  operators  on  the  tangential  vector  fields  K  and 

L  over  S. 

Under  any  orthonormal  coordinates  (u,  v)  over  S  having  w,  v  as  the  unit  basis  vectors 
and  with  n  =  u  x  v,  a  tangential  vector  field  A  over  S  can  be  written  in  matrix  form  as: 


A  = 


Then  /i  x  A  = 


-A. 


=  -io2A  where  o2  = 


0    -i 

i     0 


is  one  of  the  Pauli  spin 


matrices.  Note  that  o22  =  1.  Using  these  matrix  notations,  we  can  rewrite  Eq.  (2-12)  in  the 
following  form: 


(2-13) 


0 

io2 

K* 

-K 

M 

-N 

K 

-  2 

f?  inc 

E, 

tan 

io2 

0 

V 

-L 

N 

M 

L 

H'nc 

J 

L    J 

tan 

C.         IMPEDANCE  BOUNDARY  CONDITION 

Maxwell's  equations  alone  cannot  determine  the  electromagnetic  fields  completely. 
If  S  is  an  open  surface,  appropriate  boundary  condition  which  the  fields  satisfy  on  S  must  be 
specified.  It  is  usually  given  in  terms  of  the  impedance  boundary  condition,  a  linear  relation 

among  the  tangential  components  of  the  total  E  and  the  total  H  fields  on  S.  If  S  is  a  closed 

surface,  there  are  two  possibilities:  One  is  to  specify  the  electric  and  magnetic  properties  of 


the  volume  within  S  and  require  the  fields  to  satisfy  regularity  conditions  within  S  and  be 
linked  to  the  fields  outside  through  boundary  conditions  across  S;  another  is  to  specify  the 
impedance  boundary  condition  on  S  +  for  an  exterior  problem  or  on"  5  for  an  interior 
problem.  Note  that  an  impedance  boundary  condition  over  a  closed  surface  S  completely 
separates  the  exterior  from  the  interior  of  S.  Therefore,  the  surface  impedance  on  S~  can  be 
arbitrary  for  an  exterior  problem  while  that  on  S  +  can  be  arbitrary  for  an  interior  problem. 
In  this  thesis,  normalized  surface  impedances  Z*  are  assumed  to  be  specified  on  S  *  whether 
S  is  an  open  or  a  closed  surface.  Note  that  an  open  surface  S  can  be  considered  as  bounded 
within  the  closed  surface  formed  by  joining  S+  and  S '. 

The  impedance  boundary  conditions  on  S±  are  defined  by: 


n*  x  (£*  x  n*)  =  Z±  (n±  x  H*) 


(2-14) 


or  equivalently,  in  terms  of  the  total  surface  currents: 


«*  x  L  =  Z*  K~ 


(2-15) 


With  the  matrix  notations  for  tangential  vector  fields  over  S  in  the  orthonormal  coordinate 
system  (u,  v)  introduced  before,  we  can  consider  Z*  as  2x2  matrices  and  rewrite  Eq.  (2-15) 
in  the  form: 


+  io2  L*  =  Z*K*=  -Z±[K  ±  (K +  -  K')] 


(2-16) 


which  can  readily  be  recast  into  a  relation  among  sum  and  difference  currents: 


-A    -io. 


K*  -  K 
V  -L 


Z         0 
-A    -ia. 


(2-17) 
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with 


Z  =  -  (Z+  +  z-) 

2 


(2-18) 


and 


A   =  -  (Z  + 
2 


z-) 


(2-19) 


Eqs.  (2-13)  and  (2-17)  are  a  set  of  four  two-dimensional  vector  equations  to  be  solved  for  the 
sum  and  difference  equivalent  electric  and  magnetic  surface  current  densities  on  S. 


D.        ALGEBRA  OF  THE  SUM-DIFFERENCE  CURRENT  EQUATIONS 

The  existence  and  uniqueness  of  a  solution  to  either  the  exterior  or  the  interior 
problem  specified  in  terms  of  the  impedance  boundary  condition  have  been  well  established 
[3].  Here  we  want  to  investigate  how  such  a  solution  can  be  obtained  from  Eqs.  (2-13)  and 
(2-17).  Clearly  Eq.  (2-17)  defines  uniquely  the  algebraic  relationship  between  the  difference 
and  the  sum  currents  if  Z  is  invertible.  For  example,  the  difference  currents  can  be  expressed 
in  terms  of  the  sum  currents: 


K     -  K 
V  -  L 


0 

io2 

K 

=    - 

R 

-io2 

0 

L 

(2-20) 


where 


R  = 


Z-AZlA    -/AZ-'o. 


-io2Z  !A      o2Z  lo2 


(2-21) 
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An  equation  for  the  sum  currents  is  obtained  by  substituting  Eq.  (2-20)  into  Eq.  (2-13): 


r 

, 

r        -I 

Einc 

^tan 

M 

-N 

K 

K 

+  R 

=  2 

N 

M 

L 

L 

ffinc 

(2-22) 


which  can  be  solved  for  K  and  L.  Eq.  (2-20)  in  turn  enables  us  to  compute  the  difference 

currents  algebraically  then  split  the  sum  and  the  difference  currents  into  total  currents  on  5*. 
If  Z  is  not  invertible,  then  the  situation  is  more  complicated.  Z  can  either  be  of  rank 
0  when  Z  =  0  or  rank  1  when  det  Z  =  0  but  Z  *  0.  Eqs.  (2-13)  and  (2-17)  can  be  combined 
into  an  equation  for  the  sum  currents  K  and  L: 


1 

iA  o0 

M    -N 

K 

+ 

0 

-iZo0 

N    M 

L 

Z         0 
-A     -io. 


=  2 


1     iAo. 


0     -iZo, 


inc 
tan 

inc 

tan 


(2-23) 


When  Z  =  0,  Eq.  (2-17)  gives  the  null  relations  V  -  L~  =  io2  A  (K*  -  K~)  and 
L  =  i  o?  A  K  hence  L*  =  ;'o,  AA7*.  Eq.  (2-23)  becomes  one  for  K  only: 


[1    iAo2] 


M    -N 
N     M 


io7A 


K  =  2 


[1    /Ao2] 


inc 

tan 

inc 

tan 


(2-24) 


Eq.  (2-13)  has  to  be  used  to  find  the  difference  electric  current: 


K     -  K~  =  io2\N  +  zMo2A    K  -  2io2H™ 


(2-25) 


Therefore, 
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K+-  =  ^{\±  ia2[N  +  /Afo2A]  }£  *   io2HZC  (2-26) 


Since  the  last  term  in  Eq.  (2-26)  is  K~~     , 


£ i5C  =  i{l  ±  ia2[N  +  iMo2A]  }£  (2-27) 


L +  and  L  can  be  obtained  algebraically  by  multiplying  io2  A  to  a  and  £  respectively. 
On  the  other  hand,  by  Eq.  (2-13), 

£**  =  i*'o2{A  t  [M  -  ,Wo2A]  }tf  (2-28) 

Note  that  the  Z  =  0  case  includes  the  special  situation  Z*  =  Z'=Z  =  A  =  0  when 
both  sides  of  S  are  perfectly  conducting.  Under  this  circumstance  L  =  L~  =  0  and  the 

operator  N  is  never  involved. 

When  Z  is  rank  1 ,  Z  *  0  but  det  Z  =  0.  The  right  hand  side  of  Eq.  (2-17)  provides  one 
linear  relation  between  the  components  of  L  -  io2  AK  which  can  be  used  to  reduce  the  four 

components  of  K  and  L  as  the  unknown  quantities  in  Eq.  (2-23)  to  three  so  that  the 
remaining  three  components  can  be  solved.  The  left  hand  side  of  Eq.  (2-17)  assures  that  the 
same    linear   relationship   between    the   components    of    L  -  io2A.K    exists  _between 

corresponding  components  of  the  difference  currents.  Eq.  (2-13)  again  has  to  be  invoked  to 
compute  three  other  linearly  independent  combinations  of  the  components  of  the  difference 
currents  from  the  sum  currents. 
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E.         CONSIDERATIONS  FOR  A  CLOSED  SURFACE 

When  S  is  a  closed  surface,  the  choice  of  Z"  can  be  arbitrary  for  an  exterior  problem 
such  as  scattering  while  the  choice  of  Z+  is  arbitrary  for  an  interior  problem.  It  is  desirable 

to  choose Z"  =  -Z+  so  that  Z  =  0  and  A  =  Z+  =  -Z\  Then  we  have  L  =  io^AK  and  only  K 
has  to  be  computed.  With  such  a  choice,  for  an  exterior  problem, 


K+'sc  =  1(1  -  o2Mo2Z+  +  io2N)K  (2-29) 


£♦,«     =    i^2(Z+     +    iNo^Z+     _    M)f  (2.30) 


and  for  an  interior  problem: 


/T-"  =  1(1  -  o2Mo2Z-  -  io2N)K  (2-31) 


.«  Z°2 


L  'JL  =  — i(M  -  Z"  +  iNo0Z~)K  (2-32) 
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III.  A  THEOREM  OF  ANISOTROPIC  ABSORBERS 


A.         AXIAL    RADIATION    FROM    A    SURFACE    OF    90°    ROTATIONAL 
SYMMETRY 

The  xy-plane  cross  section  of  a  surface  5"  having  a  90  "rotational  symmetry  around 
the  z-axis  is  shown  in  Figure  3- 1 .  Because  of  this  symmetry,  S  can  be  decomposed  into  four 
non-overlapping  congruent  pieces  S,,  S2,  S3,  S4  so  that  a  90  "rotation  will  bring  5",  into  Si+]. 


Figure  3-1.  Cross  Section  of  a  Surface  of  90-Degree  Rotational  Symmetry. 


(These  subscripts  are  considered  as  equal  under  modulo  4.)  Therefore  each  piece  S.  of  S  can 
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dr.  dr. 

be  parametrized  in  the  same  orthonormal  coordinates  (u,  v),  with  u  -  — - ,  v  =  — -  the 

du  dv 

orthonormal  basis  vectors  on  Sv  as  follows: 

F,  =  (xo(u,v),yo(u,v),z0(u,v)) 

7 2  =  (-yo(u>v)>xo(u>v),zo(u,v)) 

(3-D 

r3  =  (-*0(">v)>  -y(>(u,v),z0(u,v)) 

rA  =  (yo(u,v),-xo(u,v),z0(u,v)) 

where  r.  e  Si .  As  an  example  of  a  possible  choice,  u  =  constant  and  v  =  constant  can  be  the 
lines  of  curvature  of  Sj. 

In  terms  of  the  coordinates  (u,  v),  the  sum  surface  current  distributions  on  5,  are: 


K{r)  =  Kf(u,v) 
L(f)  =  L.(m,v) 


(3-2) 


Since  for  r  »  r  , 


'*  I  ?-7„\ 

G(r)  n  £_ (3-3) 

kr 


VG(F)  -  //crG  -   -V0G(f)  ~  (3-4) 


the  radiation  from  such  current  distributions  at  a  distance  r  »  max  |  rj  |    along  the  positive 
z-axis  is,  from  Eq.  (2-6): 
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E  s\f)  =  EsJr,S,K,L) 
Jk 

4 


£-JfLpK^v)+Liy(u,v)]  (3"5) 


4 


+  yE[^v)-y«.v)] 
1=1 


ik\J(z-zo)2+x20+yl 


dudv 


Note  that  this  approximation  is  independent  of  the  wavelength;  it  is  applicable  in  regions 

closer  to  S  than  the  usual  Fresnel  zone. 

B.        CONDITION  FOR  VANISHING  ON-AXIS  BACKSCATTERING 

Consider  two  situations  when  a  linearly  polarized  plane  wave  of  unit  strength  is 
incident  on  S  along  the  z-axis  from  the  positive  direction:  Situation  1,  identified  with  the 
superscript  (1)  has  the  wave  polarized  in  the  x-direction  while  Situation  2,  identified  with  the 
superscript  (2),  has  the  wave  polarized  in  the  y-direction.  The  incident  fields  are  respectively: 

gincXD    _  £  e-ikz 

ff-.U)    =    -Se-l*  (3"6) 


gincQ.)    _    ~  e-ikz 

flinc,(2)         «      -ikz  (3-7) 

H  =  x  e    'Kz 


Note  that,  as  seen  from  the  positive  z-axis,    the  incident  wave  in  Situation  2  is  that  of 
Situation  1  rotated  by  90°  counterclockwise.  Furthermore,  Situation  2  can  be  obtained  from 

Situation  1  through  the  duality  transformation  Emc  -  H  mc,   Hmc  —  E  mc .  Therefore,  for 

a  plane  wave,  a  90°  rotation  is  equivalent  to  undergoing  the  duality  transform. 

Because  of  the  rotational  symmetry  of  S,  the  currents  excited  on  S1  under  Situation 
1  must  appear  on  Si+I  under  Situation  2.  Therefore: 
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(3-8) 


*£W) 


(3-9) 


l£\,z(u,v) 


Lu(u,v) 


Assume  that  Z  on  S  is  invertible,  the  sum  currents  on  S  are  determined  by  Eq.  (2-22).  The 
tangential  components  of  the  incident  fields  which  appear  on  the  right-hand-side  of  that 
equation  under  these  two  situations  are: 


ii-x 

ffinc,(l) 

tan 

v-x 

£r"ic,(l) 
tan 

-u-y 

-v-y 

-ikz 


(3-10) 


f?  inc,  (2) 

tan 

fjinc,(2) 
^tan 


u-y 

v-y 
u-x 
v-x 


-ikz    _ 


0     -/ 
/     0 


pine,  (I) 
^tan 

Hxzn 


(3-11) 


where  /  is  the  2x2  identity  matrix.  Therefore, 


M 

-N 

'kw 

+  R 

'kw 

=  2 

jZ  inc,  ( 1 ) 
tan 

N 

M 

P" 

L (1) 

Tjinc,(l) 
tan 

(3-12) 


M 

-N 

'gm 

+  R 

'k{2) 

=  2 

fA  inc.  (2) 
^tan 

=  2 

0 

-/ 

r?mc,(l) 
^tan 

N 

M 

L{2) 

L(2) 

frinc,(2) 

/ 

0 

ft  inc,  {I) 

1       J 

L       J 

tan 

tan 

(3-13) 


Since 


0     -/ 
/     0 

M    -N 
N     M 

- 

M    -N 
N     M 

0  -/ 

1  0 

(3-14) 


it  follows  that  if 


0     -/ 
/      0 


R  =  R 


0     -/ 
/     0 


(3-15) 


we  can  multiply 


0     -/ 
/      0 


to  Eq.  (3-12)  to  get: 


M 

-N 

0 

-/ 

KW 

+  R 

0 

-/ 

£(!)' 

=  2 

0 

-/ 

rtinc,(\) 
•^tan 

N 

M 

/ 

0 

L(1)_ 

/ 

0 

L(1)  _ 

/ 

0 

£}inc,(l) 

tan 

(3-16) 


Therefore  the  excited  surface  currents  in  these  two  situations  are  related  by: 


~K{2)' 
p2\ 

= 

0     -/ 
/     0 

= 

'-Lw 
£<i) 

(3-17) 
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Combining  this  result  with  Eqs.  (3-8)  and  (3-9),  we  have: 


K, 


K. 


a 


c- 


(2 


(2: 


:(u,v)  =  -  L™hx{utv)  =  -  K™(u,v) 

,(«,v)  =  -  L™hy(u,v)  =  K™(u,v) 

:(u,v)  =  K*{x(u,v)  =  -  L^(«,v) 

,(M,V)  =  K*\(u,v)  =  L?x\u,v) 


(3-18) 


so  that 


(i) 


£    [K?x\u,v)  +L>;)(M,v)]  =  0 


i  =  i 


(3-19) 


and 


£    [^(M,v)  -L^fov)]  =  0 


i  =  i 


Hence,  along  the  positive  z-axis,  from  Eq  (3-5), 
ik 


Es\r)-^-n\xY,[KlXXu,v)+L^u 
4iir  J  J  s\     ,=i 


,v)] 


i=i 


A^W)  -  l£\u,v) 


(3-20) 


..  r .2    2    2  (3-21) 


=  0 


and  the  backscattering  from  5"  along  the  positive  z-direction  must  vanish  if  Eq.  (3-15)  is 
satisfied. 
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C.         IMPEDANCE  MATRICES  FOR  ZERO  ON-AXIS  BACKSCATTERING 


It  can  be  verified  that  the  matrix 


Z-AZ_1A 

-iAZ_1o2 

commutes  with 

0    -/ 

-io2  Z"'A 

o0  Z'lo2 

/    0 

if  and  only  if: 


o2Z_1A  =  -AZ"'o2 


(3-22) 


Z-o2  Z~xo2   =  AZ_1A 


(3-23) 


where  both  Z  and  A  are  2  x  2  matrices.  Under  the  assumption  that  the  inverse  of  Z  exists, 
we  analyze  Eqs.  (3-22)  and  (3-23)  as  follows: 
Because  of  the  identity: 

1 


detZ 


o2ZTo2 


(3-24) 


and,  by  multiplying  o2  to  both  sides  of  Eq.  (3-22): 


Z_1  A   =  -o2  A  Z"1  o2 


(3-25) 


Eq.  (3-23)  can  now  be  transformed  to 


Z  =  -  A  (  o2  A  Z"1  o2  )  +  o2Z~l  o2 
=  -  A  o2  A  o2  (  o2  Z"1  o2)  +  o2  Z~l  o2 


/a  \2i  r?T 


detZ 


[1   -  (Ao,)2]Z 


(3-26) 
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where 


(A  o?)2  =   A  a. A  a. 


2"  ~2 


A11A22 


0 


12 


0 


A.,A22-A21 


(3-27) 


It  is  observed  that  Eq.  (3-27)  is  greatly  simplified  if  A  is  symmetric.  Therefore,  in 
this  thesis,  we  consider  only  the  case  when  A  =  Ar  .  Then   A12  =  A21  so  that: 


(Ao,)2  =  (detA)/ 


(3-28) 


From  Eq.  (3-26), 


Z  = 


1   -  det  A  ]   7T 


detZ 


(3-29) 


Eq.  (3-29)  can  be  satisfied  only  if 
Case  I 


(  1   -  detA^l 
detZ 


=  ±  1 .  These  two  cases  are: 


Z  =  -  ZT  = 


0 

-z 


12 


c12 

0 


.  Zio  *  0 


•12 


(3-30) 


detA   =  1  +  detZ  =  1  +  z 


12 


(3-31) 


CaseH 


Z  =  Z 


(3-32) 
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detZ  +  detA   =  1  (3-33) 

On  the  other  hand,  substituting  Eq.  (3-24)  into  Eq.  (3-22)  yields: 

*ll(A12-A2l)    =    ^12^2l)All  (3-34) 

z22(A12-A21)  =  (z12-z21)A22  (3-35) 

znA22  +  z22Au  =  2z21A12  =  2z12A21  (j.36) 

For  Case  1,  Eqs  (3-34)  and  (3-35)  require  that  An  =  A22  =  0,  Eq.  (3-36)  requires 
that  A,2  =  A21  =  0.  Therefore  A  =  0.  From  Eq  (3-31),  1  +  z12  =  0.  Therefore  z12  =  +  i 

so  that  Z+  =  Z~  =  Z  =  ±  o2. 

For  Case  n,  Eqs.  (3-34)  and  (3-35)  are  trivially  satisfied.  Eq  (3-36)  becomes: 

^11A22    +   *22A11    "    2*12A12    =   °  (3-37) 


Eq.  (3-33)  is  explicitly: 


^22    "    4    +    A11A22    -    An    =    1  "(3-38) 


The  sum  of  Eq.  (3-37)  with  Eq.(3-38)  is: 

(zn  +  An)(z22  +  A22)  -  (z12  +  A12)2  =  det(Z  +  A)  =  detZ+  =  1        (3-39) 
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Subtracting  Eq.  (3-37)  from  Eq.  (3-38): 

(z„  "  An)(z22  -  A22)  -  (z12  -  A12)2  =  det(Z  -  A)  =  detZ"  =  1        (3-40) 

In  summary,  two  sufficient  conditions  to  satisfy  Eqs.  (3-22)  and  (3-23)  have  been 
deduced  under  the  assumptions  that  A  is  symmetric  and  Z  is  invertible.  The  first  one  is: 

Z+  =  Z~  =  ±  o2  (3-41) 

The  other,  with  both  Z+  and  Z'  symmetric,  is: 

detZ+  =  detZ"  =  1  (3-42) 

It  should  be  noted  that  if  S  is  a  closed  surface,  the  impedance  boundary  condition 
closes  off  the  interior  of  the  surface  from  its  exterior.  Therefore  for  the  exterior  problem, 
det  Z +  =  1  is  sufficient  to  eliminate  the  on-axis  backscattering  from  a  body  of  90°  rotational 
symmetry.  Furthermore,  Z+  may  vary  with  location.  This  is  an  extension  of  Weston's  theory 
of  isotropic  absorbers  [4]. 
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IV.  SCATTERING  OF  AN  ANISOTROPICALLY  COATED  CYLINDER 

In  this  chapter,  the  sum-difference  surface  current  formulation  developed  in  Chapter 
II  is  applied  to  the  problem  of  scattering  of  an  anisotropically  coated  tubular  circular  cylinder 
of  finite  length  and  negligible  wall  thickness.  Due  to  the  rotational  symmetry,  a  Fourier  series 
expansion  can  be  utilized  to  reduce  the  variables  of  the  problem  to  only  the  one  along  the 
axis  of  symmetry,  chosen  as  the  z-axis.  The  Fourier  components  Mn ,  Nn  of  the  operators  M 

and  N  are  deduced  in  terms  of  the  Fourier  components  of  G(  r  -  r0)  and  its  partial  derivatives. 
To  solve  the  set  of  integrodifferential  equations  so  obtained,  the  equations  and  the  surface 
currents  are  both  weighted  and  expanded  over  the  Chebyshev  polynomials.  The  resultant 
infinite  system  of  linear  equations  are  then  truncated  and  inverted  numerically. 
A.         GEOMETRY  AND  COORDINATE  SYSTEM 

Figure  4- 1  shows  the  geometry  and  coordinate  system  of  a  tubular  circular  cylinder 
of  radius  a  and  length  2/z.  The  cylindrical  coordinate  system  (p,  $,  z)  is  scaled  so  that  the 
surface  of  the  cylinder  S  is  specified  by  p  =  1  and  -  1  <  z  <  1 .  The  radial  vector  is  thus  given 
by: 

f  =  app  +  hzz  (4l) 

-  a  p  (cos(J)jc  +  sin^j)  +  hzz 

To  facilitate  representing  the  tangent  vectors  over  S  in  matrix  form,  we  chose  u  =  <p  and 
v  =  z  as  the  orthonormal  basis  vectors  on  S,  so  that  n  =  p  is  the  outward  normal  to  S. 

To  properly  classify  the  polarization  of  the  incident  wave,  the  rectangular  coordinate 
system  (x,y,z)  is  determined  as  follows:  When  a  plane  wave  is  not  incident  along  the  z-axis, 
the  coordinates  are  chosen  so  that  the  wave  is  propagating  in  the  zx-plane  incident  from  the 
half-plane  containing  the  positive  x-axis.  Axial  incidences  then  are  considered  as  the  limiting 
cases  as  the  direction  of  incidence  approaches  the  positive  or  the  negative  z-axis  from  this 
half-plane.  Therefore,  even  for  axial  incidence,  a  linearly  polarized  incident  wave  having  its 
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Figure  4-1.  The  Geometry  and  Coordinate  System. 

electric  field  intensity  vector  E  pointing  in  the  y-direction  is  considered  a  TE  wave  while  one 
having  its  E  vector  in  the  zx-plane  is  considered  a  TM  wave.  In  the  spherical  coordinate 
system  (/-,0,4»),  the  direction  of  propagation  of  the  incident  wave  k  is  given  by 
k  =  -  £  cos  0 .  -  x sin  0 .  where  the  incident  angle  0 .  varies  over  the  range  0  <  0 .  <  7i .  For 

an  incident  plane  wave  of  unit  strength,  the  fields  are: 
TE  -  polariation: 


E       =  y  e' 


Hm:  =  k  x  E'n'  =  (fcos0(.  -  £sin0;.)  e 


ik'7 


TM  -  polarization: 
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■rime         *      ik»f 

H       =  y  elK 


gmc    _   fiinc   x   £   _    _^CQS0_    _    fsin0.)   e> 


where  k  =  kk.  Note  that  k»  f  =  l^cos  0 .  +  l2  p  sin  6. cos  (j>  where  /,  =  kh,  and  l^  =  ka. 
On  the  surfaces  5,  the  tangential  components  of  TE-polarized  plane  wave  are: 

E™  =  coscj)  e'k'r 

H£c  =  -  cos6;.sin<j)e'£,r~  (4-1) 

#;nc  =  -  sine,.*'**' 

and  the  tangential  components  of  TM-polarized  plane  wave  are: 

E?c  =  cos0(sin(})e'':#r 

Elznc  =  sin  6V*-'  (4-2) 

H^°  =  cos(|>e,'*,? 

B.        SCATTERED  FIELDS 

From  Eq.  (2-6),  the  scattered  fields  from  the  cylinder  are  given  in  terms  of  the  sum 
equivalent  electric  current  K  and  the  sum  magnetic  current  L: 


lxl2  k         J-iJo 

+  i   r  f271  K(ro)G(f-fo)d4>0dz0  (4-3) 

J  — w  0 
£2        J-lJo 


and 
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lxl2  k         J  Jo  °  °      °    ° 


!  r2u   r 

+ 


1  /-i/o  *  ^(r°}  G{r~ro)d^od^o  (4-4) 

AV/7n27ir^)'V^(F"F>(1)^ 
Z^2       J- Jo 


Their  components  in  cylindrical  coordinates  are: 

1       8       rl     /*2it 


/2p  3(J)  J-i  Jo        z    °  o>    v0     0 

+  lh  C  ( P  *  ^  K^o)  G{"~° #0 dZ"  (4-5) 

+  —  -^-  P  [2*  KJr)G(r-r)d<$>  dz 

,2  aP84>  J-Jo     *  °  °    *°    ° 

/,/2  dpdzJ-i  Jo       z    °  °     *°    ° 


£-  Z^(r)  =  -  f  f  / '    r  ( P  •  4>„)  L,(F  )  GCr-r^)  #o  <fe0 
/,/2  /,    oz  J-i  Jo 

/2  aP  J->  Jo      zV "  "      °    ° 

+  l  ]\  C  ( ^  W  G(F_Fo)  d*°  dZ°  (4-6) 

+  —  —  V   [2nKJr)G(f-r)d<b  dz 

/,2p  a^J-iJo     *  °         °     °   • 

+  — —  f1    (2*  K  (r)G(r-r)d<$>  dz 

lJ2p  d^dzJ-i  Jo        z    °  °        °     ° 
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i    a 


—  £/c(r)  =  -  — -^  p    r2,t($.$)L.(r)G(r-r)J(J)  <fc 
/,/2     z  /2p  dp  J -i  Jo  °     *    °  "        °      ° 

+  —-4"    P    f2n(f>^)  LJr)G(r-r)d<$>  dz 
l2p  d<$>  J-i  Jo  °      *    °  "        °     ° 


i       d2 


+  —  -?-rf     {     K.(r)G(r-r)d$  dz 

l\CKS7o)G(T-fo)d$0dz0 


1        1    #\ 

r2    3,2 


/,'  az- 


(4-7) 


iZL#;c(F)  =  -i  If    f2,I($.$)K.(r)G(r-F)^cj>  <fe 

/,/2    p  /,  dz  J-i  Jo         °    *  °         °     °    ° 

+  —  —   f    f2%  K(r)G(f-r)d§  dz 
/2p  acj>  J-i  Jo       z   °  °    *°    ° 

+  !'     '|2,(p4jVr)G(r-r)40^ 

■J  ~  I  ^  0 

+  — -^-  f1  f2,tL.(r  )G(r-F)^(J)  <fe 
+  -i~-L,f1    f2*L(r)G(r-r)dcJ>  dz 

/,/2  apazJ-i  Jo     z  °         °     °    ° 


(4-8) 


/,/2   *         /,  azJ-iio         °    *  °         °     °    ° 

l2  dp   J-i  Jo        z    "  °'    *°     ° 


/22P  a* 


i      a2 


/j/2p  5(j)5z 


f\rt<t<VG<r-V<i*o<ko 

J  - 1    J  o 


(4-9) 
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An 
1 


i     d 


i     a2 


Vi 


+  { 


-2—  f1  [2*LAr)G(r-r)d$  dz 
6z3<J)  J-iJo       *    "  °       °     ° 

J  -  1    J  0 


/,2  **, 


(4-10) 


Note  that  in  the  above  equations, 


p    =   pocos((J)-(i)o)   +  $osin((J)-4)o) 


(4-11) 


4>   =   -p0sin((J)-(})0)   +  $0cos(4>-4g 


(4-12) 


Because  of  the  rotational  symmetry  G(r-r()  depends  on  (J)  -<j)o  and  Fourier  series 

can  be  introduced  to  eliminate  the  variable  ({) .  Define  the  Fourier  series  expansion  of  a 
function/((j),z)  by: 


/«M)  =    £  ein*Uz) 


(4-13) 


then 


and 


k  \r-r       „=« 


(4-14) 
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J  -n  In 


(4-15) 


Eqs.  (4-5)  to  (4-10)  become: 
2  „«:,       N  1      a   r\ 


tKn^z)  -~tfzLL^[G^+G-^ 


-  JZL    f1  L  (z)Gdz 
/2p   J-i    zn   °     "    ° 


/22ap 


i        3    m 
/j/2    6pJ-i 


£W 


GA 


(4-16) 


/Z£*-(p'z)  =  ^  I/',  W~.-<w«.*  IA  j;  w% 


'1*2 


+<7',^ 

0) 

±<*-. 

n2 
+  G    .)— —  G 

/2p 

n         rl 

V2p  •'-' 

a 
az 

0 

G  is 

dz. 


(4-17) 


_2_ 

v2 


7I£-(P'Z) 


1       r\ 


f     hn(O[(n-l)G^-(n  +  l)Gn+l}dz0 

J  -1 


2/2p   J-i 

—  —   f  £*(z)[G    ,+G    ,]iz 
2/2  ap  J -i     ♦"   °      "~1       "  *      ° 

n     a     r  l 


GA  +'7.'WGA 


/2az  J-i 


^<o 


(4-18) 
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T^OP.:)    -±  I/.',  WG.-,  ♦<?..,]*. 


2/,    dzJ-i 


.2p 

1     n 


+  I  /_,  L*<0  [Gn_x-Gn^dz0-^j-p_x  L,n(zo)  G„<fe0 


(4-19) 


i       d    r\ 


&w 


G  dz 


i  a 


__£_  r1  [i.L  (z)]  Gdz 


(4-20) 


■JjH^z)    =-  ^L-  £  V0)[(«-DG„,-(n  +  l)GnJ^ 


lr2 

i  a 


+  _L_^_   f1  [_£.£  (Z)]c<fe    +i  [*  L  {z)  G  dz 
i2  dz  J -i    dz    z  °        J -i     z 

ii  o 


(4-21) 


Note  that  in  Eqs.  (4- 1 6)  to  (4-2 1 ),  Gn  stands  for  Gn(/j  |  z  ~z0  \ ,  /2,  P )  •  In  shifting  the  z-derivative 

to  the  current  densities  K  (z0)  and  L  (z0) ,  the  fact  that  edge  conditions  [5]  require  that  these 
components  of  the  scattering  currents  vanish  as  |  z0|  approach  1  from  below  is  used. 
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As  p  -  1*.  the  operators  Mn  and  Nn  can  be  deduced  from  Eqs.  (4-16)  to  (4-21): 


Mn.n  =  "  il\h  \\  <k0 

2      B 

P=i 

M  „  =  n  f1  dz    [G]n   .— 

M  ,.=  n—  fldz    [G]n   , 
•           8zJ-i              p 

M*,22  =    "    'V2    /'*0 

[GJp-, 

fdz1   nip=ldz0 

(4-22) 


TV 


'•21=   2/-1/Zo{[(n"1)G"-I"(n  +  1)G"  +  l]P  =  1"i/2l 


-^^-rG«-i+GnJp=i 


*U  =   ° 


(4-23) 


Note  that: 


I  dp 


-lP-r+4z\P-r\GMz-z0\,l2,p) 


dp 


JfGn(li\z-zo\>l2>1)    (4-24) 


Assuming  that  Z  is  invertible,  from  Eq.  (2-23),  the  equations  for  the  sum  currents  are 
therefore:  - 


(4-25) 


M 

n 

-n' 

n 

K 

n 

+  R 

K 

n 

=  2 

E 

tan,n 

Nn 

K 

K 

A 

tan,n 

where 
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R  = 


Z-AZ_1A     -iAZ''o. 
-/a2Z_1A      o2Z'xo^ 


(4-26) 


mc         r   "ic 


and  the  vectors Kn ,  Ln,  E^  n ,  L^    are  two  dimensional  column  matrix  representations  of 
the  respective  tangential  vector  fields  on  S  over  the  orthonomal  basis  <f),  £ ■  The  incident 
fields  are,  for  TE  and  TM  polarizations  respectively: 
TE: 


'4>,t 


n-lr'si    o;„Q\  ^  -''iZCOs6; 


H. 


n(-i)ncosQ. 


<t>." 


l2  sin  6; 


7/7         •       Q  \         -''.ZCOS0. 

7  (Z0  sin6)e     ' 


n-l  7  /-/   o;„ft  n     ~'7izcos9. 


//;;;  =  -sins.-c-irv^sine,.)* 


(4-27) 


TM: 


'<$>,n 


n(-i)ncosQt  _tt       e 

-___—  J(/  sin6.)  * 
L  sino 

2  i 

sin8/(-Onin(/2sm0f)e',7lZCOsei 


•«V»-1  jl (i    „;„ft  \  „  -i^zcosQ, 


Kn  =  (-0"-v;(/2sine,)e 


(4-28) 


C.         TRANSFORM  TO  SYSTEM  OF  LINEAR  EQUATIONS 

Since      the      surface      current      components       K.($,z0)  =0(1  ~z~)~l-       and 


2.1/2 


KA$,z0)  =  0(\  -z0)m  as  |z  |  -  1' ,  representations  of  KAzo)  and-—  K  (zj  in  Chebyshev 

dz„ 


polynomials  of  the  first  kind  combined  with  the  weighting  factor  conform  to  the  proper  edge 
behavior  of  the  currents: 
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(p,nv    o' 


1 


E  *♦!„  W 

Tisinv  4=0 


Tisinv  9=o 


E  ^n  C0S<7V 


(4-29) 


W  =  A  ECsin(^  +  1)v 


ft    «  =  o 


where  7  (z„)is  the  Chebyshev  polynomials  of  the  first  kind  and  z0  =  cosv,  -1  <z0<  1 
Similarly, 


hSZ°)=    —    EL<M   C0S<?V 

v'  7isinv  4=0 

^flW=  -  EZ,J„sin(*  +  l)v 


(4-30) 


For  an  invertible  Z,  the  coefficients  in  the  above  equations  are  to  be  determined  from 
Eq.  (4-23).  To  make  use  of  the  orthogonal  property  of  the  Chebyshev  polynomials,  the  factor 
sin  v  sin  (p  +  1 )  vfor  p  >  0  is  multiplied  to  both  sides  of  Eq.  (4-23)  before  an  integration  over 
the  range  of  v  is  carried  out.  The  results  are  described  term  by  term  in  the  subsection  to 
follow.  This  procedure  creates  an  infinite  system  of  linear  equations  to  be  solved  numerically 
after  it  is  truncated  at  an  appropriate  order  determined  by  the  electrical  size  of  the  cylinder. 


1. 


Incident  Fields 


~inc,p 
ytan,n 


H 


inc,p 
tan.n 


2     p-k  dv  ,      1A 
—  sinv  sin(/?  +  l)v 

p  +  1  Jo     71 


H, 


unc 
tan,n 

inc 
tan.n 


(4-31) 
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TE: 


47  =  ~  JtS:  Jr^COsd')  JA***>  (4-32) 


TM: 


K7  -  2*~K'  ^*iftoo.ej  /B(/2sinef) 


Zj/2sin0;     p 


£«c,P  =  (_/r+Psin6i  (  yp_i(/iC0Se.)  +  /p+1(/lCos9.))  /B(/2sin6.) 
HgP  =  {_irp-i  (  Tp^^cose,.)  +  yp  +  1(/lCos0,))  /n(/2sin9.) 


(4-33) 


where  J '{ • )  is  a  derivative  with  respect  to  the  argument.  Note  that  as  6(  approaches  0  or  n , 

only  n  =  ±1  terms  are  nonzero.  Hence  only  n  =  ±1  currents  exist.  For  axial  incident  when 
6.  =  0: 

TE: 


rinc,p    _     (~i)P     r       n\    -    irinc'P 


(4-34) 


TM: 


.wc,p   _    K~1/         j      n  \   _       c-mcp 


rjincp    _     (     QP     j        //  \    _     u'^.f 


(4-35) 


'1 
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where 


The  R  -  Matrix  Term 


— -fn —  sinv  sin(p+l)vK6n (cosv) 

P  +  \jq    n  ^ 

= HK£n  l"dv  sin(p  +  l)v  cos^v 

(/7  +  l)7I2    q-0  J0 

00 

q=0 


(4-36) 


AP'Q  =  S 


0 

P 

+q 

odd 

4 

P 

+q 

even 

7TZ(p+^  +  l)  (p-^  +  1) 


(4-37) 


/ —  sinv  sin(p  +  l)v  K (cosv) 

p  +  \Jo    7t 

;  S^n  rdv  sinv  sin(p  +  l)v  sin(<?  +  l)v  /4_38) 

0  +  1)ti2  q-o     '    Jo 

q=0 


where 


A"  =  < 


0 
-8(^  +  1) 


7r2(p+9  +  l)(/7-9  +  l)(p+9+3)(p-^-l) 


/?+<?    odd 
/?  +9    even 


(4-39) 


Therefore: 
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—  I    —  sinv  sin(/?+l)v  R 


A£q  0      0      0 

K 

n 

L 

n 

-t   * 

q   =0 

0      A™  0      0 
0      0      A£q  0 

=  £  Rp-q 

q=0 


0      0      0      Apz'q 


Kq 

n 

Lq 


K 


(4-40) 


where  Rp,q  is  a  four-by-four  matrix,  given  in  terms  of  R  by 


Ruq  ~  Ru  AS>'q 
R™  =  Ri2  A['q 

Rp»q  =  «„  Ar 


.p.? 


^iV       -    ^,4  AZ  ' 


p,q 


►     for    i  =  1,  2,  3,  4 


(4-41) 


Note  that  /?»'*  =  0  if  p+q  is  odd. 

3.         The  Ma,  Nn  Operators 

Define  the  4x4  matrix  Xpnq  by: 


2     rndv  .   ,      , , 

—  J    —  smvsin(p  +  l)v 

>  +  lJo    71 


M     -N 

n  n 


N      M 

n  n 


K 

n 

L 


q=0 


Kq 

n 


(4-42) 


— -P—  sinv  sin(p  +  l)v  M  u  [K^  (z0)] 

p  +  \J0     71 

/    — [cos/?vcos(/?+2)v]  /    cos<yvo 

q=0     p  +  l      JO     71  JO      71 

CO 

-    Y^  yP'4     Vl 
~    Z^  An,ll    A<M 

q-0 


■,  2 

-(G    ,+G    ,)-  —  G 

~v    n-l         B+l'         2      " 
to 


^<j),n 
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and, 


and, 


and, 


vP.q 

Art,ll 


p+1 


*  (riP-i  _rp*1>i  u. rP'i    r'P*2'^     n    f/-P>i    rP*2'? 

1  I 


(4-43) 


—  f  * —  sinv  sin(p+l)v  Af  „  [JST*.  (z  )] 


_s,n(p  +  l)v  —  /    cosqvoGnKln 

q=0  p  +  1     JO     71  dvJO      71 

=  f>  /"■*cos(p  +  l)v/",^cos«v„G.«;. 

9=0  JO     TC  Jo      TC 

OS 

-  Y^  yP'i   vi 

~    2s  An,21    A(}),« 
9-0 


r/?,<? 


%i  =  2n  Gn 


p+Uq 


(4-44) 


rr/o"^ sinv  sin(p+1)v  ^«  [a^(z*)] 


Jv 


£lii  r^sin(p  +  l)vi-/^cos «vJG ,.,-G nJK£n 

o=0   D  +  1    JO     7E  dWO      71 


9=0  P 

Ev/7-'?     v-1 
Ak,31    A4>,n 
9=0 


Arz,31    ~ 


.,    rrp*\,q      r,p*l,q-. 


(4-45) 
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hfo^   SlnV   Sln(/?  +  1)V  Nn,U    ^>n(Z0)] 


dv 


=  J^ — ~  f* — sinvsin(p+l)vf * — -cosqv 

a=0    D  +  \     JO      TC  JO      TC 


'    Z_/  An,41    ^fcn 
q=0 


and, 


*2  a 


(n-l)G    ,-(n  +  l)G    .-  — —  (G    ,+G    .) 


Kln 


and, 


An,41 


/, 


2(p+l)  L 


■  />.<?      /-P+2,<? 


/>.<?      r-P+2<4\ 


(n  -  1)(G™  -G^*)  -(/i  +  DCG^-G^r) 


2     oL 


(4-46) 


— -P—  sinv  sin(p  +  l)v  Mn      [K    (z0)] 

p  +  \J0     71 

=  ]C — -   f*— sinvsin(p  +  l)vru— (<?  +  l)cos(4  +  l)v  G  JT* 

^=0  P  +  l     J0     TC  JO      TC 

00 

An,12   Az,n 
4=0 


Yp,q     _        Jl(q+1)  Tr,p,q*l      rp+2,q  +  U 
Arz,12    ~        1 l^n  °*  J 


(4-47) 
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J~[fon~    SlnV    Sin^+1)V   Mn,22    \.KUZ<fi 

-    -2ill     [       dv  .       .                    Jv     .         . 
=   >    S  I    — sinvsin(/?  + l)v  /    sinv  sin(g  +  l)v  G 

q=0      p  +  l         [JO     71  JO      It 

+  _7  r~ sin(p  +  l)v— -  f71— -  (^  +  l)cos(^  +  l)v  Gw 

/2    JO     71  dvJO      71 

oo 

An,22  Az,« 
q=0 


>KZn 


and, 


Arz,22 


iliU 


)_(qP-1  +  qP+2'1  +  2  _qP>4+2  _qP  +  2>4\  _  4(<?  +  1)  ^.p  +  l.q  +  l 


P  +  l 


"(4-48) 


— -f*—  sinv  sin(p  +  l)v  AT  12  [K    (z  )] 

1  +  WO     71 


^V 


=  J^ f* — sinvsin(p  +  l)v  f* — -sinv  sin(g  +  l)v 

4=0  P+l    JO     It  io      71 


9=0  P" 

oo 

An,32   Az,n 
q=0 


K. 


and, 


An,32 


_JL  rrp''^  +  QP*24 +1  _ qP+2-1  _ qP*<1+2  "I 


4<p  +  i)  a/2 


(4-49) 


Note  that,  for  n=0,  the  expressions  simplify  to: 
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yP-1 
Ao,ll 


yPA 

Ao,41 


y/M 
Ao,22 


yP,1 
Ao,32 


yP-1 
Ao,21 


/, 


hk 


2dL 


[(Gf*  -G'( ^)+-±-^-(G™  -Gr2'q)] 


1 


p  +  1 


C(7/,'?  +  qP  +  2^+1  _qP<i  +  2  _qP^2'1\  _  £(g  +  l)f-p  +  i.<? 


j  qP-I  +  qP+2i<!+2  _qP-<]  +  2  _qP*2>1\ 


4(/?+i)  a/L 

y/>><7     _    vP.9     _    y/M     _    n 
Ao,31    ~   Ao,44    ~   Ao,12    ~    U 


(4-50) 


/>•<?_ 


Also  note  that  X0    =  0  if  p+<?  is  odd 


From  the  symmetry  of 


M     -// 


n  n 


,  we  can  deduce  that: 


yP.9     _ 

A/!,13      ~ 

An,31 

yP-V     _ 
An,23    " 

An,41 

An,33    " 

A«,ll 

y/>.<7     _ 

^■n,43    " 

yp-q 

A*,21 

yP-q     _ 
A*,14    ~ 

yp,q 

An,32 

yp.q 

An,24 

y/M 
An,42 

y/>.<7 
^n,34    ' 

yp<q 

An,12 

yP<q 

Yp,q 

An,22 

=  0 


(4-51) 


Note  that  for  p+q  odd, 


yp>q 

A/!,ll 


yP<q       .    yP<q 
An,41    "   An,22 


y/M       .    yP-4 
An,32    ~   An,23 


A*,33 


y/M       .    yP.<7      .    r> 
An,14         An,44    ""    w 


(4-52) 


and  for  p+q  even, 


yP><7 
An,21 


y/M     _    y/M     _    yP<q     _    yP<<l     _    yP-<?     _    a 
An,31    ~   A/z,12    "   A/i,13    "   AM3    ~   A«,34    ~    U 


(4-53) 


The  integrodifferential  Eq.  (4-25)  is  thus  transformed  into  the  infinite  system  of  linear 
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equations  of  the  unknown  coefficients 


Kq 

n 


L 


00 

Kq 

ginc,p 

J2[X*q+RM] 

n 

=  2 

tan.n 

<?=0 

Lq 

rjincp 

n 

tan,  n 

(4-54) 


which  is  to  be  solved  numerically. 

D.         RADIATION  IN  THE  FAR  FIELD 

In  the  far  field,  we  can  write  Eq.  (4-3)  as 

471  x 
1*2 

If      r2n 


^-Esc(r)*+ifldz0  [2nd<bo[L(ro)xr)G(r-ro) 

1,1-,  J  -1         Jo 

+  ij ' i^/2*^!  K  (r0)  -ir[  K^ro)  sin  6  sin(4»  -<j>„) 

+  ^(ro)cos6]}G(r-ro) 


(4-55) 


or  equivalently, 
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As      0  ^0 ,     only     the     n  =  ±  1     terms     are     nonzero     in     this     limit,     and 
0  =  p  =  x  cos  (J)  +jsin(j),  4>  =  -  JtsincJ)  +  jcos({).  We  have  0  ±  iq>  =  (x  ±  y)e*l^\ 


G(r)         4   fa 
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+  iy 
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(4-57) 


Similarly,  as  0  —  ti  ,  only  the  n  =  ±  1  terms  are  nonzero.  0  =  -  p   =  -fcoscf)  -  y  sincp, 
<$)  =  -x  sin  (J)  +  y  cos  4> .  We  have  0  ±  i$  =  -(x  +  y)e±l^: 
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4   , 


*"<».  <*£,■,{* 
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E.         INSIDE  AND  OUTSIDE  SURFACE  CURRENTS 

From  Eq.  (2-20)  we  have: 


K  7  -  K 


n  n 


Z"'A 


iZ~lo. 


-io2(Z-AZ-]A)    -o2AZ-1o2 


K: 


(4-59) 


Therefore  the  following  matrix  equation  is  obtained: 
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(4-60) 
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V.  COMPUTATION  AND  RESULTS 

The  solution  to  the  problem  of  the  scattering  of  an  anisotropically  coated  tubular 
cylinder  of  finite  length  as  formulated  in  the  previous  chapter  has  been  coded  in  FORTRAN 
and  tested.  The  program  listings  are  included  in  the  Appendix.  Computation  has  been  carried 
out  on  the  32-bit  Sun  SPARC  Station  running  under  the  Unix  operating  system  in  the 
Electrical  and  Computer  Engineering  Department  and  the  Computer  Center.  The  evaluation 
of  the  double  series  expansion  coefficients  of  the  Green's  function  and  its  derivatives  for 
greater  values  of  ka  and  kh  have  also  been  done  on  the  64-bit  Cray  Y-MP  EL98  in  the 
Visualization  Lab  so  that  the  accuracy  of  the  results  can  be  accessed. 

The  program  accepts  the  geometry  of  the  cylinder,  the  surface  impedance  matrices, 
the  incident  wave  frequency,  direction  and  its  polarization  in  terms  of  TE,  TM  or  some  linear 
combination  of  the  two.  It  computes  the  sum  surface  currents  and  the  far  field  radiation  in 
the  directions  specified,  including  the  strength  and  the  phase  of  each  field  component.  It  also 
breaks  down  the  sum  surface  currents  into  the  outside  and  the  inside  currents. 

In  this  chapter,  some  interesting  results  of  computation  for  the  scattering  of  a  tubular 
cylinder  having  the  length-to-diameter  ratio  hla  of  4  or  6  are  presented.  All  figures  are 
attached  at  the  end  of  this  chapter.  The  wave  is  incident  from  the  positive  z-axis  and  is 
polarized  in  the  y  direction. 
A.         COMPARISON  WITH  EXPERIMENTAL  DATA 

For  backscattering  from  a  perfectly  conducting  tubular  cylinder  of  a  y-polarized  plane 
wave  incident  from  the  positive  z-axis,  experimental  data  are  available  [6]  over  a  frequency 
band  of  well  beyond  three  octaves,  with  the  circumference-to-wavelength  ratio 
ka  =  2iza/X  varied  from  0.9448  to  3.3152.  These  data  are  measured  with  two  sets  of  four 
cylinders  each;  one  set  having  hla  =  4,  the  other  with  hla  =  6.  Both  data  sets  use  the  inner 
radii  of  the  tubular  cylinders  as  the  parameter  a  [7].  The  cylinder  length-to-wavelength  ratio 
2hlk  varies  from  1 .2030  to  4.2210  for  the  hla  =  4  case  and  from  1 .8045  to  6.3315  for  the  hla 
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=  6  case. 

The  experimental  data  are  plotted  against  the  results  of  theoretical  computation. 
Figure  5-1  shows  the  backscattering  from  the  set  of  cylinders  having  hla  =  4.  Figure  5-2 
shows  those  data  from  the  set  of  cylinders  with  hla  =  6.  The  output  of  theoretical 
computation  follows  the  measured  data  very  closely.  Note  that  the  cutoff  frequency  of  the 
dominant  circular  waveguide  mode,  TE,,,  occurs  at  2hl\  =  2.344  for  the  cylinders  with  hla 
=  4  and  at  2hl A  =  3.516  for  those  with  hla  =  6. 
B.         NULL  ON-AXIS  BACKSCATTERING 

Three  different  ways  of  coating  a  surface  impedance  Zs  on  a  tubular  cylinder  of  hla 
=  6  are  considered:  Both  on  the  outside  surface  and  the  inside  surface;  Only  on  the  inside  and 
leaving  the  outside  surface  perfectly  conducting;  Only  on  the  outside  and  leaving  the  inside 
surface  perfectly  conducting.  Here  Zs  has  the  elements  zsl]  =  0.5,  zs22  varies  from  0.1  to  5,  zsl2 
=  zs2i  =  0.  At  a  fixed  frequency  for  which  2hl  X  =  3.194,  slightly  below  the  T^  circular 
waveguide  dominant  mode  cutoff  of  3.516,  the  scattered  fields  are  plotted  in  Figures  5-3 
through  5-5.  Figure  5-3  shows  the  results  of  computation  for  the  case  Z+  =  77  =  Zs  =  Z.  As 
22  is  varied  through  2,  the  backscattering  cross  section  vanishes  as  predicted  in  Chapter  3. 
Figure  5-4  shows  the  results  for  the  case  Z+  =  0  and  Z"  =  Zs.  As  the  impedance  on  the  inside 
surface  is  increased,  the  excited  field  inside  the  tubular  cylinder  is  dissipated  and  the 
backscattered  power  decreases  exponentially.  Figure  5-5  shows  the  results  for  the  case  Z+  = 
Zs  and  Z"  =  0.  The  backscattered  power  drops  off  rapidly  at  first  as  the  impedance  on  the 
outside  surface  is  increased.  But  the  cross  section  quickly  settles  down  to  a  fixed  value 
presumably  due  to  the  current  excited  on  the  perfectly  conducting  inside  surface  of  the 
cylinder. 

Results  of  computation  for  the  same  configurations  but  at  a  higher  frequency  for 
which  2hlX  =  4.865,  above  the  TEU  circular  waveguide  dominant  mode  cutoff  of  3.516,  are 
plotted  in  Figures  5-6  through  5-8.  Now  that  the  incident  wave  can  propagate  through  the 
cylinder  in  the  dominant  waveguide  mode,  the  backscattering  cross  sections  are  about  an 
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order  of  magnitude  smaller  than  in  previous  cases.  Figure  5-6  shows  the  results  of 
computation  for  the  case  Z*  =  Z~  =Z%  =  Z.  Again  the  backscattering  cross  section  vanishes 
as  z22  is  varied  through  2.  Figure  5-7  shows  the  results  for  the  case  Z+  =  0  and  Z"  =  Zs  while 
Figure  5-8  shows  the  results  for  the  case  Z+  =  Zs  and  Z  =  0.  It  appears  that,  above  cutoff,  the 
contribution  to  the  backscattering  cross  section  from  the  inside  of  the  tubular  cylinder  is 
minimal:  once  the  outside  current  is  reduced  by  the  increase  in  surface  impedance,  the 
backscattering  cross  section  is  reduce  by  more  than  10  dB  as  shown  in  Figure  5-8.  When  the 
impedance  coating  is  applied  in  the  inside  surface,  the  maximal  reduction  in  the 
backscattering  cross  section  is  only  about  1.2  dB. 
C.        FREQUENCY  DEPENDENCE 

The  axial  backscattering  of  the  two  cases  when  the  tubular  cylinder  is  coated  only  on 
the  inside  or  only  on  the  outside  with  zslI  =  0.5  and  z^2  =  2  are  investigated  for  different 
frequencies  with  2hlX  varying  from  0.1  to  7.5.  Figure  5-9  shows  the  results  of  the  case  when 
only  the  inside  surface  is  coated  so  that  the  backscattering  is  mainly  due  to  the  current 
excited  on  the  outside  surface.  The  reflection  from  the  ends  of  the  cylinder  causes  the 
fluctuation  in  backscattering  cross  section.  Being  waves  in  free  space  on  the  outside  of  the 
cylinder,  the  maxima  and  minima  are  evenly  spaced  with  the  minima  occurring  when  2hl  X 
is  a  multiple  of  half  integer.  Figure  5-10  shows  the  results  when  only  the  outside  surface  is 
coated  and  the  current  on  the  inside  surface  dominates  the  contribution.  The  distinct  feature 
in  this  case  is  that  the  backscattering  cross  section  does  not  fluctuate  with  varying  frequency 
below  the  waveguide  mode  cutoff.  The  incident  wave  is  able  to  penetrate  deeper  into  the 
cylinder  with  increasing  frequency,  resulting  in  a  constantly  rising  strength  of  the 
backscattered  field.  Once  beyond  the  circular  waveguide  mode  cutoff,  the  wave  can  pass 
through  the  cylinder  in  the  TEn  mode  and  the  backscattering  diminishes.  The  oscillation  in 
the  cross  section  at  these  higher  frequencies  represents  the  interference  of  reflected  waves 
at  the  ends  of  the  tube  and  the  separation  between  maxima  and  minima  should  be  determined 
by  the  guide  wavelength  at  the  particular  frequency.  These  two  situations  should  be 
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compared  to  Figure  5-11  which  shows  the  results  when  both  sides  of  the  cylinder  are 
perfectly  conducting  and  the  current  can  flow  freely.  The  distinct  notch  in  the  cross  section 
near  the  TEn  mode  cutoff  at  2hl  X  =3.516  and  the  subsequent  faster  variation  in  the  cross 
section  shows  the  combination  of  the  two  distinct  features  of  Figures  5-9  and  5-10. 
D.         COMPUTATION  ACCURACY 

The  main  difficulty  encountered  in  the  computation  is  the  evaluation  of  Gp'q(lvl2) 
and  its  /2-derivative  by  double  power  series  sum  when  /,  becomes  large.  Computations  for 
Figure  5-11,5-13  and  5-15  use  the  Gpnq  values  evaluated  with  the  Cray  computer  which  has 
a  128-bit  double  precision  number.  Computations  for  Figures  5-12,  5-14  and  5-16  use  the 
Gpnq  values  evaluated  with  the  Sun  SPARC  Station  which  has  a  64-bit  double  precision 
number.  For  2hlX  greater  than  about  6.2,  the  SPARC  Station  fails  to  provide  accurate 
results. 

Figures  5- 1 3  and  5- 1 5  are  axial  backscattering  from  a  cylinder  of  hi  a  =  6  coated  with 
the  impedances  having  the  elements  zn+  =  z2{  =  0.5,  z22+  =  zu~  =  0.4,  z12+  =  z21+  =  z12"  =  z2\ 
=  0.3.  Figure  5-13  shows  the  co-polarized  backscattered  field  while  Figure  5-15  shows  the 
cross-polarized  backscattering. 
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Figure  5-6. 
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Figure  5-9. 
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Y-component  of  axial  backscattering  (perfectly  conducting  cylinder  ,  hi  a  =  6 ,  Cray  G) 
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Figure  5-11 
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Y-component  of  axial  backscattering  (perfectly  conducting  cylinder,  hi  a  =  6 ) 
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Figure  5-12. 
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Figure  5-14. 
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Figure  5-15. 
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VI.  CONCLUSIONS 

In  this  thesis,  the  sum-difference  surface  current  formulation  is  introduced  for  solving 
electromagnetic  boundary  value  problems  when  impedances  are  specified  on  both  sides  of 
a  surface.  For  an  impedance  coated  body,  the  body  can  be  treated  as  being  a  surface 
separating  the  space  into  two  regions  of  identical  medium.  For  an  exterior  problem,  the 
impedance  normalized  to  the  medium  on  the  inside  surface,  Z",  can  be  chosen  arbitrarily; 
and  for  an  interior  problem,  that  on  the  outside  surface,  Z + ,  can  be  arbitrary.  The  choice 
when  Z~  =  -Z*  is  of  particular  interest  because  the  integrodifferential  equation  has  only  the 
sum  of  the  equivalent  electric  surface  currents  on  the  outside  and  the  inside  surfaces  as  its 
unknown  to  be  solved. 

This  formulation  preserves  the  duality  nature  of  Maxwell's  equations  and  carries  it 
over  into  the  algebraic  form  of  the  integrodifferential  operators  in  the  equations  for  the  sum 
currents.  Since  a  90°  rotation  is  equivalent  to  undergoing  a  duality  transform  for  an  incident 
plane  wave,  this  particular  symmetry  in  the  algebraic  form  of  the  operators  leads  to  the 
sufficient  conditions  that  if  Z+  =  Z"  =  ±o2,  or  if  Z+  and  Z"  are  symmetric  and 
det  Z +  =  det  Z "  =  1  ,  the  on-axis  backscattering  of  an  anisotropic  impedance  coated  scatterer 
having  a  90°  rotational  symmetry  will  be  eliminated.  Note  that  in  the  symmetric  case,  Z + 
and  Z "  may  vary  with  location.  This  is  an  extension  of  Weston's  result  [4]  for  which  the 
surface  impedance  is  isotropic. 

A  FORTRAN  program  has  been  written  which  accepts  the  geometry  of  the  cylinder, 
the  surface  impedance  matrices,  the  incident  wave  frequency,  direction  and  its  polarization 
in  terms  of  TE,  TM  or  some  linear  combination  of  the  two.  It  computes  the  sum  surface 
currents  and  the  far  field  radiation  in  the  directions  specified,  including  the  strength  and  the 
phase  of  each  field  component.  It  also  breaks  down  the  sum  surface  currents  into  the  outside 
and  the  inside  currents. 

The  results  of  computation  using  this  program  agree  with  measured  data  of 
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backscattering  from  conducting  tubular  cylinders  over  a  frequency  band  of  more  than  three 
octaves.  For  a  cylinder  coated  with  surface  impedance  matrices  satisfying  the  criteria  for  null 
on-axis  backscattering,  the  numerical  computation  also  validated  the  theoretical  assertion. 
Difficulties  have  been  encountered  about  the  computational  accuracy  in  the 
evaluation  of  the  double  series  Chebyshev  expansion  coefficients  of  the  Green's  function 
G%'q(lvl2)  and  its  /2-derivative  by  double  power  series  sum  when  the  length  of  the  cylinder 
is  large  compared  to  the  wavelength:  a  compiler  with  a  64-bit  double  precision  number  can 
only  handle  a  cylinder  having  a  length  up  to  about  6.2  wavelengths.  Further  work  to  explore 
the  feasibility  of  asymptotic  evaluation  of  these  coefficients  is  recommended. 
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APPENDIX    PROGRAM  LISTING 


A. 


INCLUDE  FILES 


REALTP . INC 

C    TYPE  STATEMENTS  FOR  REAL  AND  INTEGERS  AND  DEFINITIONS  OF  CONSTANTS. 
IMPLICIT  DOUBLE  PRECISION  (A,  B,  D-H,  O-Z) 
IMPLICIT  INTEGER* 4  (I-N) 

PARAMETER  (PI  =  3 . 141592  653 5897  93  23 84 62 64 3  3  83 27D0 , PI2  =  PI  +  PI , 
+  PISQ=PI*PI) 

(ONE=l . DO , TWO=2 . DO , THR=3 . DO , HXD=16 . DO , ZERO=0 . DO , 
DEGPI=180. DO, EPS8=2. 22 044 60492503 13D-16) 
(ONEN=-ONE,HALF=ONE/TWO,THIR=ONE/THR,QUAR=HALF*HALF) 


PARAMETER 


PARAMETER 


CMPXPT . INC 

C    IMPLICIT  TYPE  STATEMENT  AND  CONSTANTS  FOR  COMPLEX  NUMBERS. 
IMPLICIT  DOUBLE  COMPLEX  (C) 

PARAMETER  ( CZERO= ( 0 . dO , 0 . dO ) , CONE= ( 1 . dO , 0 . dO ) ) 
PARAMETER  (CONEN= (-1 . dO , 0 . dO) , CI1= (0 . dO, 1 . dO )  ,  CI2= (0  .  dO  , -1 .  dO  ) 


B. 


INPUT  DATA  FILES 


CLYGEOM.PRM 
30 
5 
8 


Maximum  kh  (integer) 
Maximum  ka  (integer) 
NREGNS  (integer) 


CYLFLPT . PRM 
1024 
64 


floating  point  zero  I0BIT  (1024  bits) 
floating  point  precision  IFPBIT  (64  bits) 


INPUTDAT . PRM 
.1D-1 
.875D-2 
40 
479 
6.  DO 

1 

0 
0.D0 

6 

3 
0.D0 

1.D0 

0.DO 

12. DO 
0 


DKA0,  minimum  ka  (REAL)  in  the  computation 

DINCA,  increment  of  ka  (REAL) 

NSTR,  start  point  (INTEGER)  in  the  computation 

NEND,  end  point  (INTEGER)  in  the  computation 

RHA,  ratio  of  h  and  a-  (REAL) 

IE,  if  incident  wave  is  TE-polarized  it  is  1,  otherwise  0 

IM,  if  incident  wave  is  TM-polarized  it  is  1,  otherwise  0 

THETAI,  incident  angle  (REAL)  is  limited  to  0  to  90  degree 

NTHTAO,  total  number  (INTEGER)  of  angle  of  THATA 

NPHI,  total  number  (INTEGER)  of  angle  of  PHI0  to  be  computed 
THETAO,  initial  theta  angle  (REAL) 

DELTHO,  increment  of  theta  angle  (REAL) 

PHIA,  initial  theta  angle  (REAL) 

DELPHI,  increment  of  phi  angle  (REAL) 
IK,  set  1  to  compute  scattering  currents  on  outer  and  inner  surface 


IMPEDNCE . PRM 
0 

( .5D0, 0.D0) 
( .4D0,0.D0) 
( .3D0,0.D0) 
( .3D0,0.D0) 
( .3D0,0.D0) 


IZ,  if  perfect  conducting,  IZ=0;  otherwise,  IZ=1. 

Impedance  Z+  (phi  phi) 

Impedance  Z-  (phi  phi) 

Impedance  Z+  (phi  z) 

Impedance  Z-  (phi  z) 

Impedance  Z+  (z  phi) 
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(.3D0.0.D0)      Impedance  Z-  (z  phi] 
(.4D0,0.D0)      Impedance  Z+  (z  z) 
(.5D0,0.D0)      Impedance  Z-  (z  z) 


C.        SETUP  PROGRAM  AND  CREATED  FILES 

PROGRAM  SETUP 
P*********************************************************************** 

C  NOTES: The  format  statement  1001  need  to  be  revised  if  other  than 

C  double  precision  real  numbers  are  used  for  DKHMAX  and  DKAMAX. 

C  Statement  1001  needs  to  be  revised  if  KHMAX  or  KAMAX  exceeds 

C  3  digits. 

C  The  format  statement  1002  need  to  be  revised  if  other  than 

C  double  precision  real  numbers  are  used  for  FZERO  AND  PRECSN. 

C  Statement  1002  needs  to  be  revised  if  I0BIT  exceeds  6  digits 

C       or  IFPBIT  exceeds  4  digits. 

P*  ***********  *  ******  *  *************************************************  *  * 

C 

INCLUDE  ' REALTP . INC ' 
INCLUDE  ' CMPXTP . INC ' 
OPEN  (20,FILE='CYLGEOM.PRM' , IOSTAT=IOS, STATUS= 'OLD' ) 

IF1IOS  .NE.  0)  THEN 
WRITE(*, 9000) 

9000  FORMAT  ('Cannot  find  the  file  CYLGEOM.PRM  containing  the  ',/, 
+   'maximum  values  for  ka  and  kh,  and  the  parameter  NREGNS . ') 

STOP 

END  IF 
READ  (20,*)  KHMAX 
READ  (20,*)  KAMAX 
READ  (20,*)  NREGNS 
CLOSE  (20) 
OPEN  (20,FILE='CYLFLPT.PRM' , IOSTAT=IOS, STATUS= 'OLD' ) 

IF(IOS  .NE.  0)  THEN 
WRITE(*, 9001) 

9001  FORMAT  ('Cannot  find  the  file  CYLFLPT . PRM  containing  floating',/, 
+   '  point  zero  bit  I0BIT  and  precision  IFPBIT.') 

STOP 

END  IF 
READ (20,*)  I0BIT 
READ (20,*)  IFPBIT 
CLOSE  (20) 

OPEN  (21, FILE= ' LIMITS . INC ' , STATUS= ' UNKNOWN ' ) 
WRITE  (21,1001)  KHMAX,  KAMAX 
WRITE  (21,1002)  I0BIT,  IFPBIT,  IFPBIT 
CLOSE  (21) 

1001  FORMAT  ( 6 X, 'PARAMETER  (DKHMAX=  ',13, '.DO,  DKAMAX=  ',13, '.DO)') 

1002  FORMAT  ( 6 X, ' PARAMETER  (FZERO= ' , 16 , ' . DO ,  PRECSN= ' , 14 , ' . DO , ' , 
+        '  IFPBIT=' ,14, ' ) ' ) 

C 

C   Part  2 

DKHMAX = KHMAX 

DKAMAX=KAMAX 

ONEDEL=ONE-EPS8 

KQDIM=INT( (DKHMAX/PI) * NREGNS +ONEDEL) 

KNDIM=INT( (DKAMAX/TWO) *NREGNS+ONEDEL) 

KQDIM=MAX ( KQDIM , 1 ) 

KNDIM=MAX(KNDIM, 1) 

OPEN  (21,FILE= 'MAINDM.INC , STATUS= 'UNKNOWN' ) 

WRITE  (21,1003)  NREGNS,  KNDIM,  KQDIM 

WRITE  (21,1004) 

WRITE  (21,1005) 

WRIE  (21,1006) 

CLOSE  (21) 

1003  FORMAT  ( 6X, ' PARAMETER  (NREGNS= ' , 12 , ' ,  KNDIM= ' , 13 , 
+         ' ,  KQDIM= ' ,13, ' ) ' ) 
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1004 

FORMAT 

1005 

FORMAT 

1006 

FORMAT 

C 

C   Part  3 

(6X, 'PARAMETER  (KNDIM1=KNDIM+1 ,  KQDIM1=KQDIM+1 ) ' ) 
(6X, 'PARAMETER  (KXCRT=2*KQDIM1 ,  KCRNT=4*KQDIM1 ) ' ) 
(6X,  'PARAMETER  (MAXNG=KNDIM1 ,  MAXPEG=KQDIM/2  +  l ,  '  , 
'  MAXPOG=KQDIMl/2) ' ) 


FZERO=I0BIT 

REFC=FZERO*LOG (TWO) -LOG (PI2  ) 
REFH=HALF* (REFC-LOG (DKHMAX) ) 
REFA=HALF* (REFC-LOG (DKAMAX) ) 


ONE 


100 


KQ2=KQDIM+2 
DMXM=KQ2 

DO  100  WHILE  ( (DMXM+HALF) 
DMXM=DMXM+ONE 

CONTINUE 
MXMREG= INT ( DMXM ) 


(LOG(DMXM/DKHMAX)+ONEN)  . LT .  REFH) 


1009 


MXMSNG=INT (TWO*DKHMAX+ONEDEL) 

MXMSNG=MAX(IFPBIT,MXMSNG,KQ2) +KNDIM+1 

OPEN  ( 2 1 , FILE= ' GPQNDM . INC ' , STATUS= ' UNKNOWN ' ) 

WRITE  (21,1009)  MXMREG,  MXMSNG 

CLOSE  (21) 

FORMAT  ( 5X, 'PARAMETER  (MXMREG= ' , 14 , ' ,  MXMSNG= 

STOP 

END 


14,  ')  ') 


GPQNDM . INC 

PARAMETER  (DKHMAX=   3  0. DO, 
PARAMETER  (FZERO=   1024. DO, 


DKAMAX=    5 . DO  ) 
PRECSN=   64. DO,  IFPBIT=   64) 


LIMITS. INC 

PARAMETER  ( DKHMAX = 
PARAMETER  (FZERO= 


3  0. DO,  DKAMAX=    5. DO) 
1024. DO,  PRECSN=   64. DO,  IFPBIT=   64! 


MAINDM . INC 

PARAMETER  (NREGNS=    8,    KNDIM=    20,    KQDIM=    77) 

PARAMETER  (KNDIM1=KNDIM+1 ,    KQDIM1=KQDIM+1 ) 

PARAMETER  (KXCRT=2 *KQDIM1 ,    KCRNT=4*KQDIM1 ) 

PARAMETER  (MAXNG=KNDIM1 ,    MAXPEG=KQDIM/2+l ,    MAXPOG=KQDIMl/2 ) 


D.         PROGRAM  MAIN 


PROGRAM  MAIN 


**********■( 


*********** 


r********i 


INCLUDE     'REALTP.INC 

INCLUDE  ' CMPXTP . INC ' 

COMMON  /GCONST/  DKH, DKA, HH, HA, HSQ , DASQ, HSQN, DASQN, HHSQ, HDASQ, 
y  RAH,RAHSQ,DHA,DAH 

/INPUT1/  DKA0,DINCA,NSTR,NEND,RHA 

/INPUT2/  IE, IM,THETAI,THESINI,THECOSI,RTHEI 

/INPUT4/  IZ, IK,IS,NYSM 


COMMON 
COMMON 
COMMON 


CALL  CHKINPUT 

OPEN  (21,FILE='rhzz41.dz' , STATUS= ' UNKNOWN ' ) 

DO  8001  IS=NSTR,NEND 

DS=IS 

DKA=DKA0+DINCA*DS 

DKH=DKA*RHA 

CALL  MAXODM 

CALL  RAPQMT 

CALL  GNPQFN 

CALL  BCKSCFL(IR,CESTH,CESPH) 
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CALL  RCSPAREA(CESTH,CESPH) 
8001   CONTINUE 
CLOSE (21) 
STOP 
END 

E.        SUBROUTINE  CHKINPUT 

SUBROUTINE  CHKINPUT 

p*************************************^ 

INCLUDE  ' REALTP . INC ' 
INCLUDE  ' LIMITS . INC ' 

COMMON  /INPUT1/  DKAO , DINCA. NSTR, NEND, RHA 
COMMON  /INPUT2/  IE, IM, THETAI , THESINI , THECOSI , RTHEI 

COMMON  /INPUT3/  RTHE, RDELT, RPHI , RDELP, NTHTAO, NPHI , THESIN, THECOS , 
+  RHPI 

COMMON  /INPUT4/  IZ, IK, IS,NYSM 

OPEN(20,FILE='INPUTDAT.PRM' , IOSTAT=IOS, STATUS= 'OLD' ) 

IF(IOS  .NE.  0)  THEN 
WRITE(*,*)  'Fail  to  open  input  file  INPUTDAT . PRM ' 
STOP 

END  IF 
p*********************************************************************** 

C   Input  kh  and  ka .  These  values  are  passed  to  other 

C   parts  of  this  program  through  the  common  block  /INPUT1/. 

READ (20,*)  DKAO 

READ (20,*)  DINCA 

READ (20,*)  NSTR 

READ (20,*)  NEND 

READ (20,*)  RHA 
p+********************************************************************** 

C   Check  against  maximum  kh  and  ka  values. 
NBTW=NEND -NSTR 
DKH0=DKA0*RHA 
DKA1=DINCA*NBTW+DKA0 
DKH1=DKA1*RHA 

IF  (DKH1  .GT.  DKHMAX)  THEN 

IF  (DKA1  .GT.  DKAMAX)  THEN 
WRITE (*,*)  'Both  kh  and  ka  values  exceed  the  maximum  allowed.' 

ELSE 
WRITE (*,*)  'The  input  kh  value  exceeds  the  maximum  value  allowed.' 

END  IF 
WRITE(*,*)  'The  execution  is  terminated.' 
CLOSE (20) 
STOP 

ELSE  IF  (DKA1  . GT .  DKAMAX) THEN 
WRITE(*,*)  'The  input  ka  value  exceeds  the  maximum  value  allowed.' 
WRITE(*,*)  'The  execution  is  terminated.' 
CLOSE (20) 
STOP 

END  IF 

IF  ( (DINCH  .LT.  ZERO)  .OR.  (DINCA  . LT .  ZERO))  THEN 
WRITE  (*,*)  'The  increment  DINCH  or  DINCA  is  less  than  zero.' 
WRITE(*,*)  'The  execution  is  terminated.' 
CLOSE(20) 
STOP 

END  IF 

C   Input  the  incident  angle  (THETAI)  and  polarization  (TE  or  TM)  of  the 

C   incident  wave.  The  incident  angle  is  limited  to  0  to  90 

C   degrees.  SIN (THETAI)  and  COS (THETAI)  are  computed  also.  These 

C  parameters  are  passed  to  other  parts  of  this  program  through  the 

C   common  block  /INPUT2/. 

READ (20,*)  IE 

READ (20,*)  IM 

READ (20,*)  THETAI 
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C   Check  the  input  value  of  incident  wave 

IF  ((IE  .NE.  1)  .AND.  (IE.  NE .  0))  THEN 
WRITE  (*,*)  'Improperly  specified  the  polarization  of  incident  ', 
+   'wave,  program  is  stoped.' 
STOP 

ELSE  IF  ( (IM  .NE.  1)  .AND.  (IM  . NE .  0))  THEN 
WRITE  (*,*)  'Improperly  specified  the  polarization  of  incident  ', 
+   'wave,  program  is  stopped. ' 
CLOSE (20) 
STOP 

END  IF 

IF  ( (THETAI  .LT.  ZERO)  .OR.  (THETAI  . GT .  90.))  THEN 
WRITE  (*,*)  'Improperly  specified  incident  angle.  ', 
+   'program  is  stopped. ' 
CLOSE (20) 
STOP 

END  IF 
C  Calculate  SIN (THETAI)  and  COS (THETAI) 
RTHEI =THETAI *  PI / DEGPI 
THESINI=SIN(RTHEI) 
THECOSI=COS (RTHEI) 
p*  *  *  *  *  ***************************  *  *  *  ******  *  *  *  *  *  *  *  *  *  ********  *  * * *  *  *  *  *  *  *  *  *  * 

C  Input  the  angles  theta  and  phi  at  which  the  scattered  fields  are 

C   to  be  computed.  They  are  specified  in  terms  of  the  initial  theta 

C   (THETAO)  and  phi  (PHIO)  angles,  their  respective  increments  DELTHO 

C   and  DELPHI,  and  the  total  numbers  of  angles  NTHTAO  and  NPHI  to  be 

C   computed.  Thus  NTHTAO  and  NPHI  must  be  integers  greater  than  1.  If 

C   either  NTHTAO=0  or  NPHI=0,  no  bistatically  scattered  fields  will  be 

C  computed.  Note  that  the  scattered  electric  field  components  are 

C   computed  for  all  the  phi -angles  at  a  fixed  theta,  before  the 

C  theta-angle  is  varied.  All  angles  are  specified  in  degrees.  Theta  is 

C   limited  to  0  to  180  while  phi  is  limited  to  0  to  3  60  degrees. 

C   SIN (THETAO)  and  COS (THETAO)  are  computed  also.   These  parameters  are 

C  passed  to  other  parts  of  this  program  through  the  common  block  /INPUT3/ 

READ (20,*)  NTHTAO 

READ (20,*)  NPHI 

IF  ((NTHTAO  .LT.  0)  .OR.  (NPHI  . LT .  0))  THEN 

WRITE(*,*)  'Improperly  specified  number  of  output  angles.  ', 
+  'Program  is  stopped. ' 

STOP 

ELSE  IF  ((NTHTAO  .EQ.  0)  .OR.  (NPHI  .EQ.  0))  THEN 

CLOSE (20) 

WRITE (*,*)  'Desired  bistatic  scattered  field  direction  has  not  ', 
+  'been  (properly)  specified,  they  will  not  be  computed.' 

NTHTAO=0 

NPHI=0 

THETAO=ZERO 

DELTHO=ZERO 

PHIO=ZERO 

DELPHI=ZERO 
ELSE 

READ (20,*)  THETAO 

READ (20,*)  DELTHO 

READ (20,*)  PHIA 

READ (20,*)  DELPHI 

END  IF 

READ (20,*)  IK 
£*********************************************************************** 

C   Check  input  values. 

IF  ( (DKH0  .LE.  ZERO)  .OR.  (DKH1  . LE .  ZERO))THEN 
WRITE(*,*)  'Invalid  kh,  program  is  stopped.' 
STOP 

END  IF 

IF  ( (DKA0  .LE.  ZERO)  .OR.   (DKA1  . LE .  ZERO))THEN 
WRITE (*,*)  'Invalid  ka,  program  is  stopped.' 
STOP 

END  IF 


73 


IF  ( (DKAO  .GT.  DKHO)  .OR.  (DKA1  .GT.  DKH1))THEN 
WRITE (*,*)  'ka/kh  >  1,  program  is  stopped.' 
STOP 

END  IF 
C   Output  angle  checking  not  required: 

IF  (NTHTAO  .EQ.  0)  GO  TO  200 
C   Chcking  output  angles: 

IF  ( (THETAO  .LT.  ZERO)  .OR.   (THETAO  . GT .  DEGPI))  THEN 
WRITE (*,*)  'The  first  output  theta-angle  lies  outside  the  0  to 
+   '180  degrees  range.  Program  is  stopped. ' 
STOP 

END  IF 
THTAIF=THETAO+ (NTHTAO-1) *DELTHO 

IF  ( (THTAIF  .LT.  ZERO)  .OR.  (THTAIF  . GT .  DEGPI))  THEN 
WRITE (*,*)  'Some  of  the  specified  output  theta-angles  lie', 
+   'outside  the  0  to  180  degrees  range.  Program  is  stopped. ' 
STOP 

END  IF 
PHIMX=TWO*DEGPI 

IF  ( (PHIO  .LT.  ZERO)  .OR.  (PHIO  . GT .  PHIMX) )  THEN 
WRITE (*,*)  'The  first  output  phi-angle  lies  outside  the  0  to  ' , 
+   '3  60  degrees  range.  Program  is  stopped. ' 
STOP 

END  IF 
PHIF=PHIO+ (NPHI-1) *DELPHI 

IF  ((PHIF  .LT.  ZERO)  .OR.  (PHIF  .GT.  TWO'DEGPU)  THEN 
WRITE (*,*)  'Some  of  the  output  phi-angles  lie  outside', 
+   'the  0  to  360  degrees  range.  Program  is  stopped.' 
STOP 

END  IF 
2  00    CONTINUE 

DPI=PI/DEGPI 

RTHE=THETAO'DPI 

RPHI=PHIO*DPI 

RDELT=DELTHO*DPI 

RDELP=DELPHI*DPI 

RHPI=90.*DPI 

THESIN=SIN(RTHE) 

THECOS=COS (RTHE) 

RETURN 

END 


F.         SUBROUTINE    MAXODM 


SUBROUTINE  MAXODM 
****************************************************************** 

INCLUDE  ' REALTP . INC ' 

INCLUDE  ' MAINDM . INC ' 

COMMON  /CRNTDM/  NMAX, MXNG, IQMAX, IQMAX1 , IQMAX2 , IXCRNT , ICRNT, MXQEG  , 
v  MXQOG 

COMMON  /GCONST/  DKH, DKA, HH , HA, HSQ , DASQ, HSQN, DASQN, HHSQ, HDASQ , 
i-  RAH,RAHSQ,DHA,DAH 

COMMON  /RTHETA/  DL1COSI , DL2SINI , DLL 

COMMON  /INPUT1/  DKAO , DINCA, NSTR, NEND, RHA 

COMMON  /INPUT2/  IE, IM, THETAI , THESINI , THECOSI , RTHEI 

COMMON  /INPUT3/  RTHE, RDELT, RPHI , RDELP, NTHTAO, NPHI , THESIN, THECOS, 
k  RHPI 

ONEDEL=ONE-EPS8 

IQMAX=INT( (DKH/PI) *NREGNS+ONEDEL) 

IQMAX=MAX ( IQMAX , 1 ) 

IQMAX1=IQMAX+1 

IQMAX2=IQMAXl+2 

IXCRNT=2*IQMAX1 

ICRNT=4*IQMAX1 

MXQEG=IQMAX/2+l 

MXQOG=IQMAXl/2 
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NMAX=INT ( (DKA/TWO) *NREGNS+ONEDEL) 

NMAX=MAX(NMAX, 1) 

MXNG=NMAX+1 
C 
C   Evaluate  SIN  and  COS  functions  to  pass  along  / 

DL1C0SI=DKH*THEC0SI 

DL2SINI=DKA*THESINI 

DLL=DKH*DL2SINI 
C 
C   Evaluate  geometrical  values  to  pass  along  /GC0NS17 . 

HH=HALF*DKH 

HA=HALF*DKA 

HSQ=DKH*DKH 

DASQ=DKA*DKA 

HSQN=-HSQ 

DASQN=-DASQ 

HHSQ=QUAR*HSQ 

HDASQ=QUAR*DASQ 

RAH=DKA/DKH 

RAHSQ=RAH*RAH 

DHA=DKH*HA 

DAH=HH*HA 
C 

RETURN 

END 

G.        SUBROUTINE    RAPQMT 

SUBROUTINE  RAPQMT 
p********************************************************************** 

INCLUDE  ' REALTP . INC ' 
INCLUDE  'CMPXTP.INC 
INCLUDE  ' MAINDM . INC ' 
DIMENSION  CO( 2,2) ,  CIR(2,2) 

DIMENSION  CZSUM(2,2) ,CZDIF(2,2) , CR1 (4,4) , CR2 (4,4) , CR3 (4,4) 
DIMENSION  CRPQ1 (KCRNT, KCRNT) , CRPQ2 ( KCRNT , KCRNT ) , 
+  CRPQ3 1 ( KCRNT , KCRNT ) , CRPQ  3  2 ( KCRNT , KCRNT ) 

COMMON  /INPUT2/  IE, IM, THETAI , THESINI , THECOSI , RTHEI 
COMMON  /INPUT4/  IZ, IK, IS,NYSM 
COMMON  /XPQTMP1/  CRPQ1 , CRPQ2 
COMMON  /XPQTMP2/  CRPQ31 , CRPQ3  2 

COMMON  /CRNTDM/  NMAX, MXNG , IQMAX, IQMAX1 , IQMAX2 , IXCRNT , ICRNT , MXQEG , 
+  MXQOG 

c 

OPEN  (20,FILE='IMPEDNCE.PRM' , IOSTAT=IOS, STATUS= ' OLD '  ) 

IF  (IOS  .NE.  0)  THEN 
WRITE (*, 9000) 
9000   FORMAT  ('Cannot  find  the  file  IMPEDNCE1 . PRM  containing  the  ',/, 
+   'impedances  of  the  inner  and  outer  surfaces. ') 
STOP 

END  IF 
C 

READ (20,*)  IZ 

IF  ((IZ  .NE.  0)  .AND.   (IZ  . NE .  1))  THEN 
WRITE  (*,*)  'Improperly  specified  the  impedance  of  cylinder  ', 
+   'surface,  program  is  stopped.' 
CLOSE (20) 
STOP 

END  IF 

IF  (IZ  .EQ.  0)  THEN 
CLOSE (20) 
RETURN 
END  IF 
C  Set  up  the  matrix  when  Z  and  Delta  are  diagonal,  or  not  diagonal  but  n  <0. 
C  Initialize  the  matrix  CRPQ1 
DO  200  1=1, KCRNT 
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DO  100  J=1,KCRNT 

CRPQKI,  J)=CZERO 

CRPQ2 (I, J)=CZERO 

CRPQ31 (I, J)=CZERO 

CRPQ32 (I, J)=CZERO 
100      CONTINUE 
200    CONTINUE 

DO  400  1=1,2 
DO  300  J=l,2 

READ  (20,*)  CO (I, J) 

READ  (20,*)  CIR(I,J) 

CZSUMd,  J)  =  (CO(I,  J)+CIR(I,  J)  )  *HALF 

CZDIFd,  J)  =  (CO(I,  J)-CIR(I,  J)  )  *HALF 
3  00      CONTINUE 
40  0    CONTINUE 

CLOSE  (20) 

CZ11=CZSUM(1,1) 


CZ12=CZSUM(1 

.2) 

CZ21=CZSUM(2 

.1) 

CZ22=CZSUM(2, 

.2) 

CRDET=CZ11*CZ22- 

CZ12* 

CZ21 

CD11=CZDIF(1, 

.1) 

CD12=CZDIF(1 

,2) 

CD21=CZDIF(2 

.1) 

CD22=CZDIF(2 

.2) 

IF  (CRDET 

.EQ. 

.  CZERO)  THEK 

[ 

WRITE  (*,9001) 

STOP 

END  IF 

C  Check  the  diagonalization 

.  of  the 

impedance  matrixes 

NSYM=1 

IF  (CZ12  . 

NE. 

CZERO)  THEN 

NSYM=0 

ELSE  IF  (CZ21 

.NE. 

CZERO) 

THEN 

NSYM=0 

ELSE  IF  (CD12 

.NE. 

CZERO) 

THEN 

NSYM=0 

ELSE  IF  (CD21  .NE.  CZERO)  THEN 

NSYM=0 
END  IF 

CRDTZ  =  CONE / CRDET 

CR1 (1, 1) =CZ11- (CD11*CD11*CZ22-CD11*CD12*CZ21-CD11*CD21*CZ12 
*  +CD12*CD21*CZ11) *CRDTZ 

CR1 (1, 2) =CZ12- (CD11*CD12*CZ22-CD12*CD12*CZ21-CD11*CD22*CZ12 
*■  +CD12*CD22*CZ11)  *CRDTZ 

CR1(2,1)=CZ21- (CD11*CD21*CZ22-CD11*CD22*CZ21-CD21*CD21*CZ12 
►  +CD21*CD22*CZ11) *CRDTZ 

CR1(2,2)=CZ22- (CD12*CD21*CZ22-CD12*CD22*CZ21-CD21*CD22*CZ12 
¥  +CD22*CD22*CZ11) *CRDTZ 

CR1 (1, 3 ) = (CD12*CZ11-CD11*CZ12) *CRDTZ 

CR1 (1,4)= (CD12*CZ21-CD11*CZ22) *CRDTZ 

CR1 (2,3)= (CD22*CZ11-CD21*CZ12) *CRDTZ 

CR1 (2,4)= (CD22*CZ21-CD21*CZ22) *CRDTZ 

CR1 (3,1)=(CD11*CZ21-CD21*CZ11) *CRDTZ 

CR1 (3,2)= (CD12*CZ21-CD22*CZ11) *CRDTZ 

CR1 (4,1)=(CD11*CZ22-CD21*CZ12)*CRDTZ 

CR1 (4,2)= (CD12*CZ22-CD22*CZ12) *CRDTZ 

CR1 (3 , 3) =CZ11*CRDTZ 

CR1 (3,4)=CZ21*CRDTZ 

CR1 (4,3) =CZ12*CRDTZ 

CR1 (4 , 4) =CZ22*CRDTZ 

DO  13  00  IQE=0,MXQEG 

IQE1=IQE+1 

IQ=2*IQE 

DQ=IQ 

DQ1=IQ+1 

IQX1=4*IQ+1 
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IQX2=IQX1+1 

IQX3=IQX2+1 

IQX4=IQX3+1 

DO  1100  IPE=0,MXQEG 

IPE1=IPE+1 

IP=2*IPE 

IPX1=4*IP+1 

IPX2=IPX1+1 

IPX3=IPX2+1 

IPX4=IPX3+1 

DP=IP 

DP1=IP+1 

DBPHI=PISQ* (DP+DQ1) * (DP1-DQ) 

BPHI=4.D0/DBPHI 

DBZ=PISQ* (DP+DQ1) * (DP1-DQ) * (DP1+DQ1+ONE) * (DP-DQ1) 

BZ=-8.D0*DQ1/DB2 

CRPQ1 ( IPX1 , IQX1 ) =CR1 (1,1) *BPHI 

CRPQ1 ( IPX2 , IQX1 ) =CR1 (2,1) *BPHI 

CRPQ1 ( IPX3 , IQX1 ) =CR1 (3,1) *BPHI 

CRPQ1 (IPX4,IQX1)=CR1 (4,1) *BPHI 

CRPQ1 (IPX1,IQX2)=CR1(1,2)*BZ 

CRPQ1 ( IPX2 , IQX2 ) =CR1 ( 2 , 2 ) *BZ 

CRPQ1 ( IPX3 , IQX2 ) =CR1 ( 3 , 2 ) *BZ 

CRPQ1 (IPX4,IQX2)=CR1(4,2)*BZ 

CRPQ1 ( IPX1 , IQX3 ) =CR1 (1,3) *BPHI 

CRPQ1 ( IPX2 , I QX3 ) =CR1 (2,3) *BPHI 

CRPQ1 ( IPX3 , IQX3 ) =CR1 (3,3) *BPHI 

CRPQ1 ( IPX4 , IQX3 ) =CR1 (4,3) *BPHI 

CRPQ1 ( IPX1 , IQX4 ) =CR1 (1,4) *BZ 

CRPQ1 ( IPX2 , IQX4 ) =CR1 (2,4) *BZ 

CRPQ1 ( IPX3 , IQX4 ) =CR1 (3,4) *BZ 

CRPQ1 ( IPX4 , IQX4 ) =CR1 ( 4 , 4 ) *BZ 
1100      CONTINUE 
13  00   CONTINUE 

DO  2300  IQO=0,MXQOG 

IQ01=IQO+l 

IQ=2*IQO+l 

IQX1=4*IQ+1 

IQX2=IQX1+1 

IQX3=IQX2+1 

IQX4=IQX3+1 

DQ=IQ 

DQ1=IQ+1 

DO  2200  IPO=0,MXQOG 

IP01=IPO+l 

IP=2*IPO+l 

IPX1=4*IP+1 

IPX2=IPX1+1 

IPX3=IPX2+1 

IPX4=IPX3+1 

DP=IP 

DP1=IP+1 

DBPHI=PISQ* (DP+DQ1 ) * (DP1-DQ) 

BPHI=4.D0/DBPHI 

DBZ=PISQ* (DP+DQ1)* (DP1-DQ) * (DP1+DQ1+ONE) *(DP-DQ1) 

BZ=-8.D0*DQ1/DBZ 

CRPQ1 (IPX1,IQX1)=CR1(1,1)*BPHI 

CRPQ1 ( IPX2 , IQX1 ) =CR1 (2,1) *BPHI 

CRPQ1 ( IPX3 , IQX1 ) =CR1 (3,1) *BPHI 

CRPQ1 ( IPX4 , IQX1 ) =CR1 (4,1) 'BPHI 

CRPQ1 ( IPX1 , IQX2 ) =CR1 ( 1 , 2 ) *BZ 

CRPQ1 ( IPX2 , IQX2 ) =CR1 ( 2 , 2 ) *BZ 

CRPQ1 ( IPX3 , IQX2 ) =CR1 ( 3 , 2 ) *BZ 

CRPQ1 (IPX4,IQX2)=CR1(4,2)*BZ 

CRPQ1 ( IPX1 , IQX3 ) =CR1 (1,3) *BPHI 

CRPQ1 ( IPX2 , IQX3 ) =CR1 (2,3) *BPHI 

CRPQ1 ( IPX3 , IQX3 ) =CR1 (3,3) *BPHI 

CRPQ1 (IPX4, IQX3) =CR1 (4,3) *BPHI 
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CRPQ1 ( IPX1 , IQX4 ) =CR1 ( 1 , 4 ) *BZ 

CRPQ1 (IPX2,IQX4)=CR1(2,4) *BZ 

CRPQ1 ( IPX3 , IQX4 ) =CR1 ( 3 , 4 ) *BZ 

CRPQ1 (IPX4, IQX4)=CR1(4,4) *BZ 
2200      CONTINUE 
23  00   CONTINUE 
C  Set  up  the  matrix  Rpq  utilized  when  Z  or  Delta  is  nor  diagonal,  and  n  <  0 

CR2 (1,2)=-CR1(1,2) 

CR2(2,1)=-CR1(2,1) 

CR2 (1,3)=-CR1(1,3) 

CR2(2,4)=-CR1(2,4) 

CR2 (3,1)=-CR1 (3,1) 

CR2(4,2)=-CR1(4,2) 

CR2 (3,4)=-CRl (3,4) 

CR2(4,3)=-CR1(4,3) 

CR2(1,1)=CR1(1,1) 

CR2(2,2)=CR1(2,2) 

CR2(1,4)=CR1(1,4) 

CR2(2,3)=CR1(2,3) 

CR2 (3,2)=CR1(3,2) 

CR2(4,1)=CR1(4,1) 

CR2(3,3)=CR1(3,3) 

CR2 (4,4)=CR1 (4,4) 

DO  3300  IQE=0,MXQEG 

IQE1=IQE+1 

IQ=2*IQE 

DQ=IQ 

DQ1=IQ+1 

IQX1=4*IQ+1 

IQX2=IQX1+1 

IQX3=IQX2+1 

IQX4=IQX3+1 

DO  3100  IPE=0,MXQEG 

IPE1=IPE+1 

IP=2*IPE 

IPX1=4*IP+1 

IPX2=IPX1+1 

IPX3=IPX2+1 

IPX4=IPX3+1 

DP=IP 

DP1=IP+1 

DBPHI=PISQ* (DP+DQ1) * (DP1-DQ) 

BPHI=4.D0/DBPHI 

DBZ=PISQ* (DP+DQ1) * (DP1-DQ) * (DP1+DQ1+ONE) * (DP-DQ1) 

BZ=-8.D0*DQ1/DBZ 

CRPQ2 ( IPX1 , IQX1 ) =CR2 (1,1) *BPHI 

CRPQ2 <IPX2,IQX1)=CR2 (2,1) *BPHI 

CRPQ2 ( IPX3 , IQX1 ) =CR2 (3,1) *BPHI 

CRPQ2 (IPX4,IQX1)=CR2 (4,1)*BPHI 

CRPQ2 ( IPX1 , IQX2 ) =CR2 ( 1 , 2 ) *BZ 

CRPQ2 (IPX2,IQX2)=CR2 (2,2)*BZ 

CRPQ2 (IPX3,IQX2)=CR2(3,2) *BZ 

CRPQ2 (IPX4, IQX2)=CR2 (4,2) *BZ 

CRPQ2 ( IPX1 , IQX3 ) =CR2 (1,3) *BPHI 

CRPQ2 (IPX2,IQX3)=CR2 (2,3)*BPHI 

CRPQ2 ( IPX3 , IQX3 ) =CR2 (3,3) *BPHI 

CRPQ2 (IPX4, IQX3)=CR2 (4,3) *BPHI 

CRPQ2 (IPX1, IQX4)=CR2 (1,4) *BZ 

CRPQ2 (IPX2, IQX4)=CR2(2,4) *BZ 

CRPQ2 ( IPX3 , IQX4 ) =CR2 (3,4) *BZ 

CRPQ2 (IPX4,IQX4)=CR2 (4,4)*BZ 
3100      CONTINUE 
3  3  00   CONTINUE 

DO  4300  IQO=0,MXQOG 

IQ01=IQO+l 

IQ=2*IQO+l 

IQX1=4*IQ+1 

IQX2=IQX1+1 
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IQX3=IQX2+1 

IQX4=IQX3+1 

DQ=IQ 

DQ1=IQ+1 

DO  4200  IPO=0,MXQOG 

IP01=IPO+l 

IP=2*IPO+l 

IPX1=4*IP+1 

IPX2=IPX1+1 

IPX3=IPX2+1 

IPX4=IPX3+1 

DP=IP 

DP1=IP+1 

DBPHI=PISQ* (DP+DQ1) * (DP1-DQ) 

BPHI=4/DBPHI 

DBZ=PISQ* (DP+DQ1 ) * (DP1-DQ) * (DP1+DQ1+ONE) * (DP-DQ1 ) 

BZ=-8.*DQ1/DBZ 

CRPQ2 (IPX1,IQX1)=CR2 (1,1) *BPHI 

CRPQ2 (IPX2,IQX1)=CR2(2,1)*BPHI 

CRPQ2 ( IPX3 , IQX1 ) =CR2 (3,1) *BPHI 

CRPQ2 ( IPX 4 , IQX1 ) =CR2 (4,1) *BPHI 

CRPQ2 (IPX1,IQX2)=CR2 (1,2)*BZ 

CRPQ2 ( IPX2 , IQX2 ) =CR2 ( 2 , 2 ) *BZ 

CRPQ2 (IPX3,IQX2)=CR2 (3,2)*BZ 

CRPQ2 (IPX4, IQX2)=CR2 (4,2) *BZ 

CRPQ2 ( IPX1 , IQX3 ) =CR2 (1,3) *BPHI 

CRPQ2 (IPX2,IQX3)=CR2 (2,3) *BPHI 

CRPQ2 (IPX3,IQX3)=CR2 (3,3) *BPHI 

CRPQ2 (IPX4,IQX3)=CR2 (4,3)*BPHI 

CRPQ2 (IPX1,IQX4)=CR2 (1,4)*BZ 

CRPQ2(IPX2,IQX4)=CR2 (2,4)*BZ 

CRPQ2 (IPX3,IQX4)=CR2 (3,4)*BZ 

CRPQ2 (IPX4,IQX4)=CR2(4,4)*BZ 
4200      CONTINUE 
43  00   CONTINUE 

C  Prepare  a  mtarix  for  computing  the  scattering  currents  on  the  outer 
C   and  the  inner  surfaces 

IF  (IK  .NE.  1)  THEN 

RETURN 
END  IF 

CR3 (1,1)=-CR1(4,1) 

CR3(1,2)=-CR1(4,2) 

CR3(1,3)=-CR1(4,3) 

CR3(1,4)=-CR1(4,4) 

CR3(3,1)=CR1(2,1) 

CR3 (3,2)=CR1(2,2) 

CR3(3,3)=CR1(2,3) 

CR3(3,4)=CR1(2,4) 

CR3(2,1)=CR1(3,1) 

CR3(2,2)=CR1(3,2) 

CR3 (2,3)=CR1(3,3) 

CR3 (2,4)=CR1(3,4) 

CR3(4,1)=-CR1(1,1) 

CR3(4,2)=-CR1(1,2) 

CR3(4,3)=-CR1(1,3) 

CR3 (4,4)=-CRl(l,4) 
C 

DO  5200  IQE=0,MXQEG 

IQE1=IQE+1 

IQ=2*IQE 

DQ=IQ 

DQ1=IQ+1 

IQX1=4*IQ+1 

10X2=10X1+1 

IQX3=IQX2+1 

IQX4=IQX3+1 

DO  5100  IPE=0,MXQEG 

IPE1=IPE+1 
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IP=2*IPE 

IPX1=4*IP+1 

IPX2=IPX1+1 

IPX3=IPX2+1 

IPX4=IPX3+1 

DP=IP 

DP1=IP+1 

DBPHI=PISQ* (DP+DQ1) * (DP1-DQ) 

BPHI=4 . /DBPHI 

DBZ=PISQ* (DP+DQ1 ) * (DP1-DQ) * (DP1+DQ1+0NE) * (DP-DQ1 ) 

BZ=-8. *DQ1/DBZ 

CRPQ31 (IPX1, IQX1) = (CONE+CR3 (1,1)) *BPHI 

CRPQ3 1 ( IPX2 , IQX1 ) =CR3 (2,1) *BPHI 

CRPQ3 1 ( IPX3 , IQX1 ) =CR3 (3,1) *BPHI 

CRPQ3 1 ( IPX4 , IQX1 ) =CR3 (4,1) *BPHI 

CRPQ3 1 ( IPX1 , IQX2 ) =CR3 (1,2) *BZ 

CRPQ31 (IPX2,IQX2)=(CONE+CR3 (2,2) ) *BZ 

CRPQ31 (IPX3,IQX2)=CR3 (3,2) *BZ 

CRPQ31 (IPX4,IQX2)=CR3 (4,2) *BZ 

CRPQ31 (IPX1, IQX3)=CR3 (1,3) *BPHI 

CRPQ31 (IPX2, IQX3)=CR3 (2,3) *BPHI 

CRPQ31(IPX3,IQX3)=(CONE+CR3 (3,3) ) *BPHI 

CRPQ31 (IPX4,IQX3)=CR3 (4,3)*BPHI 

CRPQ31 (IPX1, IQX4)=CR3 (1,4) *BZ 

CRPQ31 (IPX2,IQX4)=CR3 (2,4)*BZ 

CRPQ31 (IPX3,IQX4)=CR3 (3,4)*BZ 

CRPQ31 (IPX4, IQX4) = (CONE+CR3 (4,4)) *BZ 
C 

CRPQ32 (IPX1,IQX1)=(C0NE-CR3 (1,1) ) *BPHI 

CRPQ3  2 (IPX2,IQX1)=-CR3 (2,1) *BPHI 

CRPQ32 (IPX3, IQX1)=-CR3 (3,1) *BPHI 

CRPQ3  2 (IPX4,IQX1)=-CR3 (4,1)*BPHI 

CRPQ32 (IPX1,IQX2)=-CR3 (1,2) *BZ 

CRPQ3  2 (IPX2,IQX2)=(CONE-CR3 (2,2) ) *BZ 

CRPQ32 (IPX3, IQX2)=-CR3 (3, 2) *BZ 

CRPQ3  2 (IPX4,IQX2)=-CR3 (4,2)*BZ 

CRPQ32 (IPX1,IQX3)=-CR3 (1,3)*BPHI 

CRPQ32 (IPX2,IQX3)=-CR3 (2,3) *BPHI 

CRPQ32 (IPX3,IQX3)=(CONE-CR3 (3,3) ) *BPHI 

CRPQ32 (IPX4,IQX3)=-CR3 (4,3)*BPHI 

CRPQ3  2 (IPX1, IQX4) =-CR3 (1,4) *BZ 

CRPQ3  2 (IPX2,IQX4)=-CR3 (2,4) *BZ 

CRPQ32 (IPX3 , IQX4) =-CR3 (3,4) *BZ 

CRPQ32 (IPX4,IQX4)=(CONE-CR3 (4,4)) *BZ 
5100      CONTINUE 
52  00   CONTINUE 
C 

DO  6200  IQO=0,MXQOG 

IQ01=IQO+l 

IQ=2*IQO+l 

DQ=IQ 

DQ1=IQ+1 

IQX1=4*IQ+1 

IQX2=IQX1+1 

10X3=10X2+1 

IQX4=IQX3+1 

DO  6100  IPO=0,MXQOG 

IPEl=IPO+l 

IP=2*IPO+l 

IPX1=4*IP+1 

IPX2=IPX1+1 

IPX3=IPX2+1 

IPX4=IPX3+1 

DP=IP 

DP1=IP+1 

DBPHI=PISQ* (DP+DQ1) * (DP1-DQ) 

BPHI=4. /DBPHI 

DBZ=PISQ* (DP+DQ1) * (DP1-DQ) * (DP1+DQ1+ONE) * (DP-DQ1) 
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6100 
6200 

9001 


BZ=-8.*DQ1/DBZ 
CRPQ31 (IPXl.IQXl 
CRPQ31 (IPX2, IQX1 
CRPQ31(IPX3,IQX1 
CRPQ31 (IPX4, IQX1 
CRPQ31 (IPX1, IQX2 
CRPQ31 (IPX2, IQX2 
CRPQ31(IPX3,IQX2 
CRPQ31 (IPX4, IQX2 
CRPQ31(IPX1,IQX3 
CRPQ31 (IPX2, IQX3 
CRPQ31 (IPX3,IQX3 
CRPQ31 (IPX4,IQX3 
CRPQ31 (IPX1, IQX4 
CRPQ31 (IPX2, IQX4 
CRPQ31(IPX3,IQX4 
CRPQ31 (IPX4, IQX4 

CRPQ32 (IPXl.IQXl 
CRPQ3  2 (IPX2, IQX1 
CRPQ32 (IPX3, IQX1 
CRPQ32(IPX4,IQX1 
CRPQ32 (IPX1.IQX2 
CRPQ32 (IPX2,IQX2 
CRPQ32 (IPX3, IQX2 
CRPQ32 (IPX4,IQX2 
CRPQ32 (IPX1,IQX3 
CRPQ32 (IPX2,IQX3 
CRPQ32 (IPX3,IQX3 
CRPQ32 (IPX4, IQX3 
CRPQ32 (IPX1, IQX4 
CRPQ32 (IPX2,IQX4 
CRPQ32 (IPX3,IQX4 
CRPQ3  2(IPX4,IQX4 

CONTINUE 
CONTINUE 
RETURN 

FORMAT ( ' The  given 

►        '  sum.  Exe 

END 


= (CONE+CR3 (1,1) ) *BPHI 

=CR3 (2,1) *BPHI 

=CR3 (3,1)*BPHI 

=CR3 (4,1)*BPHI 

=CR3 (1,2) *BZ 

MCONE+CR3 (2,2) ) *BZ 

=CR3 (3,2)*BZ 

=CR3 (4,2)*BZ 

=CR3 (1,3)*BPHI 

=CR3 (2,3)*BPHI 

= (CONE+CR3 (3,3)) *BPHI 

=CR3 (4,3) *BPHI 

=CR3 (1,4) *BZ 

=CR3 (2,4) *BZ 

=CR3 (3,4) *BZ 

=(CONE+CR3 (4,4) ) *BZ 

= ( CONE-CR3 (1,1)) *BPHI 
=-CR3 (2,1) *BPHI 
=-CR3 (3,1)*BPHI 
=-CR3 (4,1)*BPHI 
=-CR3 (1,2) *BZ 
= (CONE-CR3 (2,2) ) *BZ 
=-CR3 (3,2)*BZ 
=  -CR3 (4,2) *BZ 
=  -CR3 (1,3)*BPHI 
=  -CR3 (2,3) *BPHI 
=(CONE-CR3 (3,3) ) *BPHI 
=-CR3 (4,3)*BPHI 
=-CR3 (1,4)*BZ 
=-CR3 (2,4)*BZ 
=-CR3 (3,4) *BZ 
=(CONE-CR3 (4,4) ) *BZ 


inside  and  outside  impedance  have  a  singular' 
cution  is  terminated.') 


H.        SUBROUTINE  GNPQFN 


SUBROUTINE  GNPQFN 


t *****  i 


INCLUDE  ' REALTP . INC ' 

INCLUDE  ' MAINDM . INC ' 

COMMON  /GCONST/  DKH , DKA, HH, HA, HSQ, DASQ, HSQN, DASQN, HHSQ, HDASQ, 
+  RAH,RAHSQ,DHA,DAH 

CHARACTER  FILEEVEN1*12 ,  FILEODDl*12 
C  Set  up  the  file  names 

NC1=INT(DKA) 

DC=NC1 

NC=100000 . *DKH+NC1 

DC1=DKA-DC 

NC2=INT(1000.*DC1) 
C  Set  up  the  indices  of  file  names  of  Green's  function,  and  check  whether 
C  these  files  exist.  If  these  files  exist,  then  use  it  directly.  Otherwise, 
C  call  subroutine  XPQINI1 . 

FILEEVEN1='E 

FILE0DD1='0 

IGE=8*2* (MAXPEG+1 ) * (MAXPEG+2 ) 

IGO=8*2* (MAXPOG+1) * (MAXPOG+2) 

WRITE  (FILEEVEN1 (2:8) , ' (17 .7) • )  NC 

WRITE  (FILEODD1 (2:8) ,' (17.7) ' )  NC 

WRITE  (FILEEVEN1(10:12)  ,  '  (13  .3)  '  )  NC2 
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WRITE  (FILE0DD1 (10:12)  ,'  (13.3)  '  )  NC2 

OPEN  (2  8,ACCESS='DIRECT' , FILE=FILEEVEN1 , RECL=IGE, IOSTAT=IOS, 
+  STATUS= ' OLD ' ) 

IF  (IOS  .NE.  0)  THEN 
CLOSE  (28) 
CALL  XPQINI1 (DKH,DKA) 

ELSE 
OPEN  (2  9,ACCESS='DIRECT' , FILE=FILEODDl , RECL=IGO, IOSTAT=IOS, 
+  STATUS = 'OLD' ) 

IF  (IOS  .NE.  0)  THEN 
CLOSE  (29) 
CALL  XPQINI1(DKH,DKA) 

ELSE 
CLOSE  (28) 
CLOSE  (29) 

CALL  XPQINI (DKH.DKA) 
END  IF 
ENDIF 
RETURN 
END 
C 

SUBROUTINE  XPQINI (DKHIN, DRAIN) 

INCLUDE  ' REALTP . INC ' 
INCLUDE  ' CMPXTP . INC ' 
INCLUDE  ' MAINDM . INC ' 

DIMENSION  GNE(4, (MAXPEG+1) * (MAXPEG+2) 12)  , 
+  GNO(4, (MAXPOG+1) * (MAXPOG+2) 12) 

DIMENSION  CGNE ( 0 : MAXPEG , 0 : MAXPEG, KNDIM1+1 ) , 
+  CGNO(0:MAXPOG,0:MAXPOG,KNDIM1+1) 

DIMENSION  CDGNE ( 0 : MAXPEG , 0 : MAXPEG, KNDIM1+1 ) , 
+  CDGNO (0:MAXPOG,0:MAXPOG,KNDIM1+1) 

COMMON  /CRNTDM/  NMAX, MXNG, IQMAX , IQMAX1 , IQMAX2 , IXCRNT , ICRNT , MXQEG  , 
+  MXQOG 

COMMON  /GPQTMP/  CGNE , CDGNE , CGNO , CDGNO 
SAVE  /GPQTMP/ 
C 

CHARACTER  FILEEVEN1*12 ,  FILEODDl*12 
DKH=DKHIN 
DKA=DKAIN 
c 

C  Initialize  the  matrix  CGNE,  CGNO,  CDGNE,  CDGNO 
DO  105  IC=  1,KNDIM1+1 

DO  102  IB=  0, MAXPEG 
DO  101  IA=  0,MXPEG 
CGNE ( IA, IB , IC) =CZERO 
CDGNE (IA, IB, IC)=CZERO 

101  CONTINUE 

102  CONTINUE 

DO  104  ID=0,MAXPOG 

DO  103  IE=0,MAXPOG 
CGNO (IE, ID, IC)=CZERO 
CDGNO (IE, ID, IC) =CZERO 

103  CONTINUE 

104  CONTINUE 

105  CONTINUE 
C 

C  Set  up  the  file  name 
NC1=INT(DKA) 
DC=NC1 

NC=100000 . *DKH+NC1 
DC1=DKA-DC 
NC2=INT(1000. *DC1) 
FILEEVEN1='E 
FILEODDl='0 

IGE=8*2* (MAXPEG+1) * (MAXPEG+2) 
IGO=8*2* (MAXPOG+1) * (MAXPOG+2) 
WRITE  (FILEEVEN1 (2:8) , ' (17.7) ' )  NC 
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WRITE  (FILE0DD1(2:8)  ,  '  (17.7)  '  )  NC 

WRITE  (FILEEVEN1 (10:12)  ,'  (13.3)  '  )  NC2 

WRITE  (FILEODD1(10:12)  ,  '  (13.3)  '  )  NC2 

OPEN  (2  8,ACCESS= 'DIRECT' , FILE=FILEEVEN1 , RECL=IGE, IOSTAT=IOS, 
+  STATUS='OLD' ) 

C 

OPEN  (2  9,ACCESS='DIRECT' , FILE=FILEODDl , RECL=IGO, IOSTAT=IOS, 
+  STATUS= ' OLD ' ) 

DO  900  NI=1,MXNG+1 
C     The  following  values  N,  DN  and  DNH  are  passed  to  the  G-computation 
C      related  subroutines  through  the  common  block  /NCONST/ . 

READ  (2  8,REC=NI)  GNE 

IRECE=0 

DO  500  IQE=0,MXQEG 

DO  300  IPE=0,IQE-1 

CGNE ( IPE, IQE, NI)=CGNE( IQE, IPE,NI) 

CDGNE ( IPE , IQE , NI ) =CDGNE ( IQE , IPE , NI ) 
3  00  CONTINUE 

DO  400  IPE=IQE,MXQEG 

IRECE=IRECE+1 

GR=GNE(1,IRECE) 

GI=GNE(2,IRECE) 

GDR=GNE(3, IRECE) 

GDI=GNE(4,IRECE) 

CGNE ( IPE , IQE , NI ) =DCMPLX (GR , GI ) 

CDGNE ( I PE , IQE , NI ) =DCMPLX ( GDR , GDI ) 
400  CONTINUE 

500       CONTINUE 

READ  (2  9,REC=NI)  GNO 

IRECO=0 

DO  800  IQO=0,MXQOG 

DO  600  IPO=0,IQO-1 

CGNO (IPO, IQO, NI)=CGNO( IQO, IPO, NI) 

CDGNO ( I PO , IQO , NI ) =CDGNO ( IQO , IPO , NI ) 
600  CONTINUE 

DO  700  IPO=IQO,MXQOG 

IRECO=IRECO+l 

GR=GNO(l,IRECO) 

GI=GNO(2,IRECO) 

GDR=GNO(3,IRECO) 

GDI=GNO(4, IRECO) 

CGNO ( IPO , IQO , NI ) =DCMPLX (GR ,  GI ) 

CDGNO ( I PO , IQO , NI ) =DCMPLX ( GDR , GDI ) 
700  CONTINUE 

800       CONTINUE 
900    CONTINUE 

CLOSE  (28) 

CLOSE  (29) 

RETURN 

END 

C 

SUBROUTINE  XPQINI1 (DKHIN, DKAIN) 

INCLUDE  ' REALTP . INC ' 

INCLUDE  ' CMPXTP . INC ' 

INCLUDE  ' MAINDM . INC ' 

INCLUDE  ' GPQNDM . INC ' 

INCLUDE  ' LIMITS . INC ' 

DIMENSION  GNE (4,  (MAXPEG+1 ) * (MAXPEG+2 ) /2 )  , 
+  GNO (4, (MAXPOG+1) * (MAXPOG+2) 17) 

DIMENSION  CGNE ( 0 : MAXPEG , 0 :MAXPEG, KNDIM1+1 ) , 
+  CGNO ( 0 : MAXPOG , 0 : MAXPOG , KNDIM1+1 ) 

DIMENSION  CDGNE ( 0 : MAXPEG , 0 : MAXPEG , KNDIM1+1 ) , 
+  CDGNO (0 : MAXPOG, 0 : MAXPOG, KNDIM1+1 ) 

COMMON  /CRNTDM/  NMAX, MXNG, IQMAX, IQMAX1 , IQMAX2 , IXCRNT , ICRNT, MXQEG , 
+  MXQOG 
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COMMON  /GCONST/  DKH, DKA, HH, HA, HSQ, DASQ, HSQN, DASQN, HHSQ, HDASQ, 
+  RAH,RAHSQ,DHA,DAH 

COMMON  /SUMLMT/  DHMAX, DAMAX, DSNG, IHMAX, IAMAX, ISNG 

COMMON  /NCONST/  DN,DNH,N 

COMMON  /GPQTMP/  CGNE , CDGNE , CGNO , CDGNO 

CHARACTER  FILEEVEN1*12 ,  FILEODDl*12 

SAVE  /GPQTMP/, /GCONST/, /SUMLMT/ 
C 

DKH=DKHIN 

DKA=DKAIN 
C  Initialize  the  matrix  CGNE,  CGNO,  CDGNE,  CDGNO 

DO  105  IC=  1,KNDIM1+1 

DO  102  IB=  0,MAXPEG 
DO  101  IA=  0,MXPEG 

CGNE ( IA, IB , IC ) =CZERO 

CDGNE ( IA , IB , IC ) =CZERO 

101  CONTINUE 

102  CONTINUE 

DO  104  ID=0,MAXPOG 

DO  103  IE=0,MAXPOG 
CGNO (IE, ID, IC) =CZERO 
CDGNO (IE, ID, IC) =CZERO 

103  CONTINUE 

104  CONTINUE 

105  CONTINUE 
C 

C   Determine  the  maximum  number  of  terms  for  r-  and  m-  series  sums  and 

C  pass  them  through  /SUMLMT/: 

C 

REFC=FZERO*LOG(TWO)-LOG(PI2) -ONE 

REFH=HALF* (REFC-LOG (DKH) ) 

REFA=HALF* (REFC-LOG (DKA) ) 
C 

IQMX2=IQMAX+2 

DHMAX=IQMX2 

DO  1100  WHILE  ( (DHMAX+HALF) * (LOG (DHMAX /DKH) +ONEN)  . LT .  REFH) 

DHMAX=DHMAX+ONE 
1100       CONTINUE 

IHMAX=INT (DHMAX) 
C 

DAMAX=AINT (DKA) +ONE 

DO  1200  WHILE  ( (DAMAX+HALF) * (LOG (DAMAX/DKA) +ONEN)  . LT .  REFA) 

DAMAX=DAMAX+ONE 
1200       CONTINUE 
C 

IAMAX= INT ( DAMAX ) 
C 

ISNG=INT(TWO*DKH+ONE-EPS8) 

ISNG=MAX(IFPBIT, ISNG, IQMX2 ) 
C   Checking  dimensions. 

IF (IHMAX  .GT.  MXMREG)  THEN 

WRITE(*,*)  'Warning:  IHMAX  =  ', IHMAX, '  >  MXMREG  =  ', MXMREG 

WRITE (*,*)   'IHMAX  IS  SET  TO  MXMREG  IN  XPQINI1 ' 

IHMAX=MXMREG 
END  IF 

ISN1=MXMSNG-N-1 

IF (ISNG  .GT.  ISN1)  THEN 

WRITEC*)  'Warning:  ISNG  =  MSNG,'  >  MXMSNG-N-1  =  '  ,  ISN1 

WRITE!*,*)   'ISNG  IS  REDUCED  TO  MXMSNG-N-1  IN  XPQINI1 ' 

ISNG=ISN1 
END  IF 
C 

DSNG = ISNG 
C 

FILEEVEN1='E 

FILEODDl='0 

NC1=INT(DKA) 

DC=NC1 
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NC=DKH*100000 . +NC1 

DC1=DKA-DC 

NC2=INT(1000.*DC1) 

IGE=8*2* (MAXPEG+1)* (MAXPEG+2) 

IGO=8*2* (MAXPOG+1) * (MAXPOG+2) 

WRITE  (FILEEVEN1 (2:8) ,' (17.7) ' )  NC 

WRITE  (FILE0DD1 (2:8)  ,  '  (17  .7)  ' )  NC 

WRITE  (FILEEVEN1 (10:12)  ,  '  (13 .3)  ' )  NC2 

WRITE  (FILEODD1 (10:12) ,  '  (13.3)  ' )  NC2 

OPEN  (2  8,ACCESS='DIRECT' , FILE=FILEEVEN1 , RECL=IGE, STATUS= 'UNKNOWN' ) 

OPEN  (2  9,ACCESS='DIRECT' , FILE=FILEODDl , RECL=IGO, STATUS= ' UNKNOWN ' ) 
C   Initialize  CGNE,  CDGNE,  CGNO  and  CDGNO 

DO  1900  NI=1,MXNG+1 
C     The  following  values  N,  DN  and  DNH  are  passed  to  the  G-computation 
C     related  subroutines  through  the  common  block  /NCONST/ . 

N=NI-1 

DN=N 

DNH=DN+HALF 

IRECE=0 

DO  1500  IQE=0,MXQEG 

IQ=2*IQE 

DO  1300  IPE=0,IQE-1 

IP=2*IPE 

CGNE ( IPE , IQE , NI ) =CGNE ( IQE , IPE , NI ) 

CDGNE ( IPE , IQE , NI ) =CDGNE ( IQE , IPE , NI ) 
13  00         CONTINUE 

DO  14  00  IPE=IQE,MXQEG 

IRECE=IRECE+1 

IP=2*IPE 

CALL  GDN(IP,IQ,GI,GDI,GR,GDR) 

GNE(1,IRECE)=GR 

GNE(2,IRECE)=GI 

GNE(3,IRECE)=GDR 

GNE(4,IRECE)=GDI 

CGNE (IPE, IQE,NI)=DCMPLX(GR,GI) 

CDGNE (IPE,IQE,NI)=  DCMPLX ( GDR , GDI ) 
1400         CONTINUE 
1500      CONTINUE 

WRITE  (28,REC=NI)  GNE 

IRECO=0 

DO  1800  IQO=0,MXQOG 

IQ=2*IQO+l 

DO  1600  IPO=0,IQO-1 

IP=2*IPO+l 

CGNO ( IPO , IQO , NI ) =CGNO ( IQO , IPO ,  NI ) 

CDGNO ( IPO , IQO , NI ) =CDGNO ( IQO , IPO , NI ) 
1600         CONTINUE 

DO  1700  IPO=IQO,MXQOG 

IRECO=IRECO+l 

IP=2*IPO+l 

CALL  GDN(IP,IQ,GI,GDI,GR,GDR) 

GNO(l, IRECO)=GR 

GNO(2, IRECO)=GI 

GNO(3,IRECO)=GDR 

GNO(4,IRECO) =GDI 

CGNO ( IPO , IQO , NI ) =DCMPLX (GR , GI ) 

CDGNO ( I PO , IQO , NI ) =DCMPLX ( GDR , GDI ) 
17  00         CONTINUE 
1800      CONTINUE 

WRITE  (29,REC=NI)  GNO 
1900   CONTINUE 

CLOSE (28) 

CLOSE (29) 

RETURN 

END 
C 

SUBROUTINE  GDN ( IPIN, IQIN , GIOUT , GDIOUT , GROUT , GDROUT ) 
p********************************************************************** 
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C  This  subrouitne  sets  up  P  and  Q  dependent  constants  and  passes  along 
C  /PCONST/,  then  calls  subroutines  GDREG  and  GDSNG  to  compute  G(P,Q,N) 
C   and  its  derivative. 

C 

INCLUDE  ' REALTP . INC ' 

COMMON  /PCONST/  DP, DQ, DS , DD, DSH, DDH, DSSQ, DDSQ , IP, IQ, IS , ID 
C  Check  and  transform  input  variables. 

IP=IPIN 

IQ=IQIN 

ISPQ=IP+IQ 

IS=ISPQ/2 

IF    ((ISPQ+D/2     .GT.     IS)    THEN 

GIOUT=ZERO 

GDIOUT=ZERO 

GROUT=ZERO 

GDROUT=ZERO 

WRITE  (*,  9001)  ISPQ 
9001   FORMATf'The  parameters  P  and  Q  have  a  sum  of  the  ODD  integer  ',14, 
+  /,',  G(P,  Q,  N)  and  its  derivative  have  been  set  to  0.') 

RETURN 

END  IF 

IF  (IP  .LT.  IQ)  THEN 

IP=IQIN 

IQ=IPIN 
END  IF 
C 

DP=IP 

DQ=IQ 
C 

ID=IS-IQ 

DS  =  IS 

DD=ID 

DSH=DS+HALF 

DDH=DD+HALF 

DSSQ=DS*DS 

DDSQ=DD*DD 
C 

CALL  GDREG (GI, GDI, GRR,GDRR) 

CALL  GDSNG (GRS,GDRS) 

GIOUT=GI 

GDIOUT=GDI 

GROUT=GRR+GRS 

GDROUT=GDRR+GDRS 

RETURN 

END 
CC 

SUBROUTINE  GDREG (GIOUT , GDIOUT , GROUT , GDROUT ) 

C   This  subrouitne  computes  the  regular  part  of  G{P,Q,N)  and  its 
C   derivative, 
p**************************************************** 

c 

INCLUDE  ' REALTP . INC ' 
INCLUDE  'GPQNDM.INC 

COMMON  /GCONST/  DKH, DKA, HH, HA, HSQ, DASQ, HSQN, DASQN, HHSQ,  HDASQ, 
+  RAH,RAHSQ,DHA,DAH 

COMMON  /NCONST/  DN,DNH,N 

COMMON  /PCONST/  DP, DQ, DS , DD, DSH, DDH, DSSQ, DDSQ, IP , IQ,  IS  ,  ID 

COMMON  /SUMLMT/  DHMAX, DAMAX, DSNG, IHMAX, IAMAX, ISNG 

SAVE  /GCONST/, /SUMLMT/ 
C 
C   Reserve  working  space  to  store  r-independent  numbers. 

DIMENSION  GIM(MXMREG) , GRRM (MXMREG) 
C 

C   Computation  starts. 
C 
C   Compute  overal  constant  factors: 
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GRRF=QUAR*DKH/PISQ/ (DSSQ-QUAR) / (QUAR-DDSQ) / (DN+ONE) 

GIF=HALF/DSH 

DJN=ZERO 

DO  3  00  JN=1,N 
DJN=DJN+ONE 
HATN=HA/DJN 
GRRF=GRRF*HATN*HATN 
GIF=GIF*HATN*HA/ (DJN+DSH) 

3  00       CONTINUE 

DJP=ZERO 

DO  400  JQ=1, IQ 
DJP=DJP+ONE 
GIF= (HHSQ/DJP/DJP) *GIF 

4  00       CONTINUE 

DO  500  JP=IQ+1,IP 

DJP=DJP+ONE 

GIF=(HH/DJP)*GIF 
500       CONTINUE 

IF  ((ID+D/2  .GT.  ID/2)  GIF=-GIF 
C 

C   Compute  GI,  GDI,  GRREG,  and  GDRREG . 
C     Compute  r-independent  factors  and  store  in  GIM  and  GRRM: 

SMH=DSH+ONEN 

SM1=DS 

SM2=DS+DS 

PM1=DP 

QM1=DQ 

THRH=HALF 

DMl=ZERO 

THRS=HALF+DS 

THRD=HALF+DD 

THRDN=HALF-DD 

THRSN=HALF-DS 

DO  600  JM=1,IHMAX 

SMH=SMH+ONE 

SMl=SMl+ONE 

SM2=SM2+ONE 

PMl=PMl+ONE 

QMl  =  QMl+ONE 

THRH=THRH+ONE 

DMl=DMl+ONE 

THRS=THRS+ONE 

THRD=THRD+ONE 

THRDN=THRDN+ONE 

THRSN=THRSN+ONE 

GIM(JM) = (SMH/SM2) * (SMH/PM1) * (SM1/QM1) * (HSQN/DM1) 

GRRM(JM)=(THRH/THRS) * (HSQN/THRD) * (DM1/THRDN) * (DM1/THRSN) 
600       CONTINUE 
C 

C     Compute  r-  and  m-sum  for  GI  and  GRR. 
C        Setup  initial  r  related  values 

DJR=DAMAX 

DNR=DN+DJR 

DNHR=DNR-HALF 

D2NR=DN+DNR 

DNSHR=DNR+DSH 

DN2R=DNR+ONE 
C  m-sum  for  r=IAMAX 

DNSHRM=DNSHR+DHMAX 

DN2RM=DN2R+DHMAX 

ANNEXTR=ONE+GIM(IHMAX) /DNSHRM 

ANR0R=ONE+GRRM(IHMAX) /DN2RM 
DO  700  JM=IHMAX-1, 1, -1 

DNSHRM=DNSHRM+ONEN 

DN2RM=DN2RM+ONEN 

ANNEXTR=ONE+ANNEXTR*GIM(JM) /DNSHRM 

ANR0R=ONE+ANR0R*GRRM(JM) /DN2RM 
7  00       CONTINUE 
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ANNEXT=ANNEXTR 

DANNEXT = DNR  *  ANNEXTR 

ANRO=ANROR 

DANRO=DNR*ANROR 

DO  900  JR=IAMAX,1,-1 
C        Compute  factors  for  r=JR 

ATMP= (DNHR/D2NR) * (DASQN/DJR) 

ATMPI=ATMP/DNSHR 

ATMPR=ATMP/DN2R 
C        Setup  new  r  related  values  for  r=JR-l 

DJR=DJR+ONEN 

DNR=DNR+ONEN 

DNHR=DNHR+ONEN 

D2NR=D2NR+ONEN 

DNSHR=DNSHR+ONEN 

DN2R=DN2R+ONEN 
C        m-sum  for  r=JR-l 

DNSHRM=DNSHR+DHMAX 

DN2RM=DN2R+DHMAX 

ANNEXTR=ONE+GIM(IHMAX) /DNSHRM 

ANR0R=ONE+GRRM(IHMAX) /DN2RM 
DO  800  JM=IHMAX-1,1, -1 

DNSHRM=DNSHRM+ONEN 

DN2RM=DN2RM+ONEN 

ANNEXTR=ONE+ANNEXTR*GIM(JM) /DNSHRM 

ANR0R=ONE+ANR0R*GRRM(JM) /DN2RM 
800       CONTINUE 

ANNEXT=ANNEXTR+ANNEXT*ATMPI 

DANNEXT=DNR*ANNEXTR+DANNEXT*ATMPI 

ANR0=ANR0R+ANR0*ATMPR 

DANR0=DNR*ANR0R+DANR0*ATMPR 
900       CONTINUE 
C 

GIOUT=ANNEXT*GIF 

GDIOUT=DANNEXT*GIF 

GROUT=ANR0  *GRRF 

GDROUT=DANR0  *GRRF 
C 

RETURN 

END 
C 

SUBROUT INE  GDSNG ( GROUT , GDROUT ) 

C   This  subrouitne  computes  the  'singular'  part  of  G(P,Q,N)  and  its 

C   derivative, 
p***************************************************^ 

c 

INCLUDE  ' REALTP . INC ' 

INCLUDE  ' GPQNDM . INC ' 

COMMON  /GCONST/  DKH, DKA, HH , HA, HSQ , DASQ, HSQN, DASQN, HHSQ, HDASQ, 
+  RAH,RAHSQ,DHA,DAH 

COMMON  /NCONST/  DN,DNH,N 

COMMON  /PCONST/  DP, DQ, DS , DD, DSH, DDH, DSSQ, DDSQ, IP,  IQ,  IS  ,  ID 

COMMON  /SUMLMT/  DHMAX, DAMAX, DSNG , IHMAX, IAMAX, ISNG 

SAVE  /GCONST/, /SUMLMT/ 
C 
C   Reserve  working  space  to  store  r-independent  numbers. 

DIMENSION  GSIM(MXMSNG) , GS2M (MXMSNG) , 
+  GR2M (MXMSNG) , GR2R (0 : MXMSNG) , GDR2R (0 : MXMSNG) 

C 

C   Computation  starts. 
C 

MXGR2M=ISNG+N 
c 
C  Compute  overal  constant  factors: 

GRSF=QUAR/PISQ/DKH 
C 
C   Compute  GR1,  GDR1 ,  GR2 ,  and  GDR2 : 
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C     Compute  initial  constants: 

S1LHA=TW0*L0G (HXD/RAH) 

DK=HALF 

S1SDO=ZERO 

DO  1000  JK=1,ID 

S1SDO=S1SDO+ONE/DK 

DK=DK+ONE 
1000      CONTINUE 

SlSD0=TWO*SlSD0 

DO  1100  JK=ID+1,  IS 

SlSD0=SlSD0+ONE/DK 

DK=DK+ONE 
1100      CONTINUE 

SlSD0=-TWO*SlSD0 

SN0N=ZERO 

SN10=ZERO 

SN2  0=ZERO 

IF  (N  .GT.  0)  THEN 

DJN=ONE 

DJNH=HALF 

DO  1200  JN=1,N-1 

DJNI=ONE/DJN 

DJNHI=ONE/DJNH 

SN10=DJNI*DJNHI+SN10 

SN20=DJNHI*DJNHI* ( (ONE-QUAR*DJNI) *HALF*DJNI+ONE) +SN20 

DJN=DJN+ONE 

DJNH=DJNH+ONE 
1200         CONTINUE 

DJNI=ONE/DJN 

DJNHI=ONE/DJNH 

SN0N=S1LHA+HALF*DJNHI -SN1 0  *QUAR 

SN10= (DJNI*DJNHI+SN10) *QUAR 

SN20=HALF* (DJNHI*DJNHI* ( (QUAR*DJNI+ONEN) *HALF*DJNI+ONEN) -SN20) 
END  IF 

SN10=S1LHA-SN10 

SN20=PISQ*THIR+SN20 
C 
C      Compute  r-independent  terms  and  store  in  GR2M,  GS1M,  and  GS2M: 

GR2JM=ONE 

S1SDM=S1SD0 

S2SDM=ZERO 

DM=ZERO 

DMH=-HALF 

DO  13  00  JM=1,MXGR2M 

DM=DM+ONE 

DMH=DMH+ONE 

DMI=ONE/DM 

DMHI=ONE/DMH 

DMHSQ=DMH*DMH 

DMHSS=DMHSQ-DSSQ 

DMHSD=DMHSQ-DDSQ 

GR2JM=(DMHSS*DMI*DMI) * (DMHSD*DMI*DMHI ) *GR2JM 

GR2M(JM) =GR2JM 

TWSMHS=TWO*DSSQ/DMHSS 

TWDMHD=TWO*DDSQ/DMHSD 

S1SDM=S1SDM- (TWSMHS+TWDMHD+DMI+ONE) *DMHI 

S2SDM=S2SDM- ( (TWSMHS+THR) *TWSMHS+ (TWDMHD+THR) *TWDMHD+ 
+        (ONE-QUAR*DMI) *TWO*DMI+ONE) /DMHSQ 

GS1M(JM)=S1SDM 

GS2M(JM)=S2SDM 
13  00      CONTINUE 
C 

C      Compute  m-sum  for  GR1 ,  GDR1  and  store  in  GR2R,  GDR2R 
IF  (N  .GT.  0)  THEN 

DJR=DN+ONEN 

DJRH=DJR-HALF 

SN0R=SN0N 
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SGR1M=SN0R+S1SD0 

SDGR1M=DJR*SGR1M+0NEN 

DJM=ONE 

DJRM=DJR 

DO  1400  JM=1,N-1 

SR1M=SN0R+GS1M ( JM) 

HSQF=HSQ/DJM/DJRM 

SGR1M=SGR1M*HSQF+SR1M*GR2M(JM) 

SDGR1M=SDGR1M*HSQF+ (DJR*SR1M+0NEN) *GR2M(JM) 

DJM=DJM+ONE 

DJRK=DJRM+ONEN 
1400         CONTINUE 

GR2R (N-l) =SGR1M*TW0 

GDR2R (N-l ) =SDGR1M*TW0 

DO  1600  JR=N-2,0,-l 

SN0R=SN0R+ONE/DJRH+ONE/ (DN-DJR) -ONE/ (DN+DJR) 

DJR=DJR+ONEN 

DJRH=DJRH+ONEN 

SGR1M=SN0R+S1SD0 

SDGR1M=DJR*SGR1M+0NEN 

DJM=ONE 

DJRM=DJR 

DO  1500  JM=1,JR 

SR1M=SN0R+GS1M ( JM) 

HSQF=HSQ/DJM/DJRM 

SGR1M=SGR1M*HSQF+SR1M*GR2M ( JM) 

SDGR1M=SDGR1M*HSQF+ (DJR'SRIM+ONEN) *GR2M(JM) 

DJM=DJM+ONE 

DJRM=DJRM+ONEN 
1500  CONTINUE 

GR2R ( JR) =SGRlM*TWO 

GDR2R ( JR) =SDGRlM*TWO 
1600         CONTINUE 

END  IF 
C 
C     Compute  m-sum  for  GR2 ,  GDR2  and  store  in  GR2R,  GDR2R 

DRN=DN 

SN1R=SN10 

SN2R=SN20 

SR2MM=SN1R+S1SD0 

SR2RMS=SR2MM*SR2MM+SN2R 

DR2RMS=DRN*SR2RMS-TWO*SR2MM 
DJM=ONE 

DJRM=DRN 

DO  1700  JM=1,N 

SR2MM=SN1R+GS1M ( JM) 

SR2RM=SR2MM*SR2MM+SN2R-tG.92M(JM) 

DR2RM=DRN*SR2RM-TWO*SR2MM 

HSQF=HSQ/DJM/DJRM 

SR2RMS=SR2RMS*HSQF+SR2RM*GR2M(JM) 

DR2RMS=DR2RMS*HSQF+DR2RM*GR2M ( JM) 

DJM=DJM+ONE 

DJRM=DJRM+ONEN 
17  00   CONTINUE 

GR2R(N)=SR2RMS 

GDR2R(N)=DR2RMS 

DJR=ZERO 

DO  1900  JR=N+1,MXGR2M 

DJR=DJR+ONE 

DRN=DRN+ONE 

DJRI=ONE/DJR 

DRNHI=ONE/ (DRN-HALF) 

DR2NI=ONE/ (DRN+DN) 

DNHRF=DNH*DRNHI*DR2NI 

SN1R=DJRI-DNHRF+SN1R 

SN2R=DJRI*DJRI- (DRNHI+DR2NI ) *DNHRF+SN2R 

SR2MM=SN1R+S1SD0 

SR2RMS=SR2MM*SR2MM+SN2R 
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DR2RMS=DRN*SR2RMS-TWO*SR2MM 

DJM=ONE 

DJRM=DRN 

DO  1800  JM=1,JR 

SR2MM=SN1R+GS1M ( JM) 

SR2RM=SR2MM*SR2MM+SN2R+GS2M(JM) 

DR2RM=DRN*SR2RM-TWO*SR2MM 

HSQF=HSQ/DJM/DJRM 

SR2RMS=SR2RMS*HSQF+SR2RM*GR2M(JM) 

DR2RMS=DR2RMS*HSQF+DR2RM*GR2M ( JM) 

DJM=DJM+ONE 

DJRM=DJRM+ONEN 
1800         CONTINUE 

GR2R(JR)=SR2RMS 

GDR2R(JR)=DR2RMS 
1900      CONTINUE 
C 

C        Compute  r-sum: 
C 

IF  (RAHSQ  .GT.  HALF)  THEN 
C 

C  Convergence  acceleration  via  Euler-Abel  Transf omation : 

C 

DJR=ZERO 

DJRH=DJR-HALF 

DJRN=DJR+DN 

DJRNN=DJR-DN 

IF  (N  .GT.  0)  THEN 

RFCTOR=ONE/DN 

GR2R ( 0 ) =RFCTOR*GR2R ( 0 ) 

GDR2R ( 0 ) =RFCTOR*GDR2R ( 0 ) 
DO  2000  JR=1,N-1 

DJR=DJR+ONE 

DJRH=DJRH+ONE 

DJRN=DJRN+ONE 

DJRNN=DJRNN+ONE 

RFCTOR= (DJR*DJRH/DJRN/DJRNN) *RAHSQ*RFCTOR 

GR2R ( JR) =RFCTOR*GR2R ( JR) 

GDR2R(JR)=RFCTOR*GDR2R(JR) 
2000  CONTINUE 

DJR=DJR+ONE 

DJRH=DJRH+ONE 

DJRN=DJRN+ONE 

DJRNN=DJRNN+ONE 

RFCTOR=HALF* (HALF-DN) *RAHSQ*RFCTOR 

GR2R (N) =RFCTOR*GR2R (N) 

GDR2R (N) =RFCTOR*GDR2R (N) 
ELSE 

RFCTOR=ONE 

END  IF 

DO  2100  JR=N+1,MXGR2M 

DJR=DJR+ONE 

DJRH=DJRH+ONE 

DJRN=DJRN+ONE 

DJRNN=DJRNN+ONE 

RFCTOR= (DJR*DJRH/DJRN/DJRNN) *RAHSQ*RFCTOR 

GR2R ( JR) =RFCTOR*GR2R ( JR) 

GDR2R ( JR) =RFCTOR*GDR2R ( JR) 
2100         CONTINUE 
C 

C  Compute  zeroth-order  delta-r  term: 

DO  2300  JDELTR=1,MXGR2M 

DO  2200  JR=MXGR2M, JDELTR, -1 

GR2R ( JR) =GR2R ( JR-1 ) -GR2R ( JR) 

GDR2R(JR) =GDR2R(JR-1) -GDR2R(JR) 
2200  CONTINUE 

23  00         CONTINUE 

SGR2R=HALF*GR2R (MXGR2M) 
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SDGR2R=HALF*GDR2R (MXGR2M) 

DO    2400    JR=MXGR2M-1, 0, -1 

SGR2R= (SGR2R+GR2R(JR) ) *HALF 

SDGR2R= (SDGR2R+GDR2R(JR) ) *HALF 
240  0  CONTINUE 

GROUT=GRSF*SGR2R 

GDROUT=GRSF*SDGR2R 
C 

ELSE 
C 

C  Direct  summation  when  RAHSQ  is  no  greater  than  1/2: 

C 

DJR=DSNG 

DJRN=DJR+DN 

D JRNH = D JRN - HALF 

DJR2N=DJRN+DN 

RFCTOR= (DJRN*DJRNH/DJR/DJR2N) 'RAHSQ 

SGR2R=RFCTOR*GR2R (MXGR2M) 

SDGR2R=RFCTOR*GDR2R (MXGR2M) 

DO  2500  JR=MXGR2M-1,N+1, -1 

DJR=DJR+ONEN 

DJRN=DJRN+ONEN 

DJRNH=DJRNH+ONEN 

DJR2N=DJR2N+ONEN 

RFCTOR= (DJRN*DJRNH/DJR/DJR2N) *RAHSQ 

SGR2R= (GR2R(JR) -SGR2R) *RFCTOR 

SDGR2R= (GDR2R ( JR) -SDGR2R) *RFCTOR 
2500         CONTINUE 

IF  (N  .EQ.  0)  THEN 

GROUT=(GR2R(0) -SGR2R) *GRSF 

GDROUT=(GDR2R(0) -SDGR2R) *GRSF 
ELSE 

RFCTOR=HALF* (HALF-DN) *RAHSQ 

SGR2R= (GR2R (N) -SGR2R) *RFCTOR 

SDGR2R= (GDR2R (N) -SDGR2R) *RFCTOR 

DJR=DN 

DJRH=DJR-HALF 

DJRN=DJR+DN 

DJRNN=DJR-DN 

DO  2600  JR=N-1, 1, -1 

DJR=DJR+ONEN 

DJRH=DJRH+ONEN 

DJRN=DJRN+ONEN 

DJRNN=DJRNN+ONEN 

RFCTOR= (DJR*DJRH/DJRN/DJRNN) * RAHSQ 

SGR2R=(GR2R(JR)-SGR2R) *RFCTOR 

SDGR2R= (GDR2R ( JR) -SDGR2R) *RFCTOR 
2  600  CONTINUE 

GROUT= (GR2R ( 0 ) -SGR2R) *GRSF/DN 

GDROUT= (GDR2R ( 0 ) -SDGR2R) 'GRSF/DN 
END  IF 
END  IF 

RETURN 

END 


I.  SUBROUTINE  BCKSCFL 

SUBROUTINE  BCKSCFL ( IR , CESTH , CESPH ) 

C 

INCLUDE  ' REALTP . INC ' 

INCLUDE  ' CMPXTP . INC ' 

INCLUDE  'MAINDM.INC 
C 

REAL*  4  AKPHPR . AKPHPI , AKZPR , AKZPI , ALPHPR , ALPHPI , ALZPR ,  ALZPI , 
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+        AKPHNR , AKPHNI , AKZNR , AKZNI , ALPHNR , ALPHNI , ALZNR , ALZNI 

DIMENSION  CIW(KCRNT) ,CKLN(KCRNT) 

DIMENSION  CIWO(KXCRT) , CKLNO (KXCRT) , CXPQNO (KXCRT, KXCRT) 

DIMENSION  CRPQ1 ( KCRNT , KCRNT ) , CRPQ2 (KCRNT, KCRNT) , 
+  CXPQN ( KCRNT , KCRNT ) , CXRPQ ( KCRNT , KCRNT ) 

DIMENSION  CKPHP (61, 61) ,CKZP(61, 61) , CLPHP (61, 61) , 
+  CLZP(61,61),  CKPHN(61, 61) ,CKZN(61, 61)  , 

+  CLPHN(61,61) ,CLZN(61,61) 

COMMON  /INPUT2/  IE, IM, THETAI , THESINI , THECOSI , RTHEI 

COMMON  /INPUT4/  IZ, IK, IS,NYSM 

COMMON  /NCONST/  DN,DNH,N 

COMMON  /XPQTMP/  CXPQN, CXRPQ, CXPQNO 

COMMON  /XPQTMP1/  CRPQ1 , CRPQ2 

COMMON  /CKLMTX/  CKPHP, CKZP, CLPHP, CLZP, CKPHN, CKZN, CLPHN, CLZN 

COMMON  /CRNTDM/  NMAX, MXNG , IQMAX, IQMAX1 , IQMAX2 , IXCRNT, ICRNT, MXQEG , 
+  MXQOG 

COMMON  /INPUT3/  RTHE, RDELT, RPHI , RDELP, NTHTAO, NPHI , THESIN, THECOS , 
+  RHPI 

CESTH=CZERO 

CESPH=CZERO 

DO  100  JQ=1, KXCRT 

CKLNO (JQ)=CZERO 
10  0    CONTINUE 

DO  200  JQ=1, KCRNT 

CKLN(JQ) =CZERO 

2  00    CONTINUE 

DO  400  IF=-30,30 

DO  300  JZ=-30,30 
CKPHP ( IF , JZ ) =CZERO 
CKZP (IF, JZ)=CZERO 
CLPHP ( IF , JZ ) =CZERO 
CLZP (IF, JZ)=CZERO 
CKPHN ( IF , JZ ) =CZERO 
CKZN (IF, JZ)=CZERO 
CLPHN ( IF , JZ ) =CZERO 
CLZN (IF, JZ)=CZERO 

3  00       CONTINUE 

4  00    CONTINUE 
C 

IF  (RTHEI  .EQ.  0 . DO )  GO  TO  2000 
C  If  incident  angle  is  not  equal  to  0,  use  this  loop 

DO  1800  NA=0,NMAX 
C 

CALL  INCIDNT(NA,CIW0,CIW) 
C 

IF  (NA  .EQ.  0)  THEN 
CALL  XPQ0 

ELSE 
CALL  XPQN(NA) 
ENDIF 
C  If  the  cylinder  is  made  of  perfect  conductor  and  no  coating  on  it 
IF  (IZ  .EQ.  0)  THEN 
DN=NA 
C  Use  IMSL  library  to  solve  linear  system 

CALL  DLSACG ( IXCRNT , CXPQNO , KXCRT , CIW0 , 1 , CKLNO ) 
CALL  ESCFAR (NA, CKLNO , CKLN , CETHN , CEPHN ) 
C 

C   If  we  calculate  in  X  and  Y  components  as  theta  approaches  0  or  pi, 
C   then  theta  component  is  cosin  phi  in  X  direction  plus  sine  phi  in  Y 
C   direction,  and  phi  component  is  negative  sine  phi  in  X  direction  plus 
C   cosin  phi  in  Y  direction. 

IF  (RTHE  .EQ.  0 . DO )   THEN 
CETHN=COS(RPHI) 'CETHN+SIN (RPHI ) *CEPHN 
CEPHN=-SIN(RPHI) *CETHN+COS (RPHI ) *CEPHN 

ELSEIF  (RTHE  .EQ.  PI)  THEN 
CETHN= -COS (RPHI ) *CETHN-SIN (RPHI ) *CEPHN 
CEPHN=-SIN(RPHI) *CETHN+COS (RPHI ) *CEPHN 
END  IF 
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CESTH=CESTH+CETHN 
CESPH=CESPH+CEPHN 

IF  (NA  .NE.  0)  THEN 

IF  (IE  .EQ.  1)  THEN 
DO  600  I=2,KXCRT,2 
CKLNO (I) =-CKLN0 (I) 
60  0    CONTINUE 

END  IF 

IF  (IM  .EQ.  1)  THEN 
DO  800  I=1,KXCRT,2 
CKLN0(I)=-CKLN0(I) 
800    CONTINUE 

END  IF 
NA1=-NA 
DN=NA1 
CALL  ESCFAR (NA1 , CKLNO , CKLN, CETHN, CEPHN) 

IF  (RTHE  .EQ.  0 . DO )   THEN 
CETHN=COS (RPHI ) *CETHN+SIN (RPHI ) *CEPHN 
CEPHN=-SIN (RPHI ) *CETHN+COS (RPHI ) *CEPHN 

ELSEIF  (RTHE  . EQ .  PI)  THEN 
CETHN=- COS (RPHI) *CETHN-SIN (RPHI) 'CEPHN 
CEPHN= -SIN (RPHI) *CETHN+COS (RPHI ) * CEPHN 
END  IF 
C 

CESTH=CESTH+CETHN 
CESPH=CESPH+CEPHN 
END  IF 
END  IF 
C  If  the  cylinder  is  coated  with  anisotropic  material 
IF  (IZ  .EQ.  1)  THEN 
DO  1100  I=1,KCRNT 

DO  1000  J=1,KCRNT 
CXRPQ ( I , J) =CXPQN ( I , J) +CRPQ1 ( I , J) 
1000     CONTINUE 
1100   CONTINUE 
DN=NA 
CALL  DLSACGdCRNT,  CXRPQ,  KCRNT,CIW,  1,CKLN) 

IF  ((IK  .EQ.  1)  .AND.  (IR  . EQ .  45))  THEN 
CALL  KLCRNT(DN,CKLN,CIW) 

END  IF 
CALL  ESCFAR (NA, CKLNO , CKLN, CETHN, CEPHN) 

IF  (RTHE  .EQ.  0 . DO )   THEN 
CETHN=COS (RPHI ) *CETHN+SIN (RPHI ) *CEPHN 
CEPHN=-SIN(RPHI) *CETHN+COS (RPHI ) *CEPHN 

ELSEIF  (RTHE  . EQ .  PI)  THEN 
CETHN=- CCS (RPHI) *CETHN-SIN (RPHI ) *CEPHN 
CEPHN=-SIN(RPHI) *CETHN+COS (RPHI ) *CEPHN 
END  IF 
C 

CESTH=CESTH+CETHN 
CESPH=CESPH+CEPHN 

IF  (NA  .NE.  0)  THEN 

IF  (NYSM  .EQ.  0)  THEN 
DO  1300  I=1,KCRNT 

DO  1200  J=1,KCRNT 
CXRPQ ( I , J) =CXPQN (I , J) +CRPQ2 (I , J) 
1200      CONTINUE 
13  00   CONTINUE 

ELSE 
DO  1500  I=1,KCRNT 

DO  1400  J=1,KCRNT 
CXRPQ (I, J)=CXPQN(I, J) +CRPQ1 (I, J) 
1400     CONTINUE 
1500   CONTINUE 

END  IF 
CALL  DLSACG(ICRNT,CXRPQ,KCRNT,CIW, 1,CKLN) 
IF  (IE  .EQ.  1)  THEN 
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DO  1600  I=2,KCRNT,4 
11=1+1 

CKLN<I)=-CKLN(I) 
CKLN(I1)=-CKLN(I1) 
1600   CONTINUE 

END  IF 

IF  (IM  .EQ.  1)  THEN 
DO  1700  I=1,KCRNT,4 
12=1+3 

CKLN(I)=-CKLN(I) 
CKLN(I2)=-CKLN(I2) 
17  00   CONTINUE 

END  IF 
NA1=-NA 
DN=NA1 

IF  ((IK  .EQ.  1)  .AND.  (IS  .EQ.  199))  THEN 
CALL  KLCRNT(DN,CKLN,CIW) 

END  IF 
CALL  ESCFAR (NA1 , CKLNO , CKLN , CETHN , CEPHN) 

IF  (RTHE  .EQ.  0.D0)   THEN 
CETHN=COS(RPHI) *CETHN+SIN (RPHI ) "CEPHN 
CEPHN=-SIN(RPHI) *CETHN+COS (RPHI ) "CEPHN 

ELSEIF  (RTHE  .EQ.  PI)  THEN 
CETHN= -COS (RPHI ) *CETHN-SIN (RPHI ) 'CEPHN 
CEPHN=-SIN(RPHI) *CETHN+COS (RPHI) * CEPHN 
END  IF 
C 

CESTH=CESTH+CETHN 
CESPH=CESPH+CEPHN 
END  IF 
END  IF 
1800   CONTINUE 

C  If  the  cylinder  is  coated  with  anisotropic  material  and  we  want  to  calculate 
C  equivalent  currents  K  nad  L  on  inner  and  outer  surfaces. 
IF  (IK  .EQ.  1)  THEN 

IF  (IS  .EQ.  199)  THEN 
OPEN(31,FILE='khzz41.dz ' , STATUS =' UNKNOWN ' ) 
OPEN(32,FILE='lhzz41.dz ' , STATUS= 'UNKNOWN1 ) 

END  IF 
DO  1900  IF=-30,30 

DO  1850  JZ=-30,30 
AKPHPR=REAL(CKPHP(IF, JZ) ) 
AKPHPI=IMAG(CKPHP(IF, JZ) ) 
AKZPR=REAL(CKZP(IF, JZ) ) 
AKZPI=IMAG(CKZP(IF, JZ) ) 
AKPHNR=REAL ( CKPHN ( IF , JZ ) ) 
AKPHNI=IMAG(CKPHN(IF, JZ) ) 
AKZNR=REAL(CKZN(IF, JZ) ) 
AKZNI=IMAG(CKZN(IF, JZ) ) 
ALPHPR=REAL(CLPHP(IF, JZ) ) 
ALPHPI=IMAG(CLPHP(IF, JZ) ) 
ALZPR=REAL(CLZP(IF, JZ) ) 
ALZPI=IMAG(CLZP(IF, JZ) ) 
ALPHNR=REAL(CLPHN(IF, JZ) ) 
ALPHNI=IMAG(CLPHN(IF, JZ) ) 
ALZNR=REAL(CLZN(IF, JZ) ) 
ALZNI=IMAG(CLZN(IF, JZ) ) 


WRITE (31, *) 

AKPHPR, ' 

' ,AKPHPI 

WRITE (31, *) 

AKZPR, ' 

' ,AKZPI 

WRITE (32, *) 

ALPHPR, ' 

' ,ALPHPI 

WRITE (32, *) 

ALZPR, ' 

1 ,ALZPI 

WRITE (31, *) 

AKPHNR, ' 

' ,AKPHNI 

WRITE (31, *) 

AKZNR, ' 

' ,AKZNI 

WRITE (32, *) 

ALPHNR, ' 

' , ALPHNI 

WRITE(32, *) 

ALZNR, ' 

' ,ALZNI 

1850 

CONTINUE 

1900 

CONTINUE 
CLOSE (31) 
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CLOSE (32) 

END  IF 
RETURN 
C  If  incident  angle  is  equal  to  zero,  the  program  should  go  to  this  loop 
C  because  it  need  to  compute  when  n=+l  and  n=-l  only. 
2000   NA=1 
DN=NA 

CALL  INCIDNT(NA,CIWO,CIW) 
CALL  XPQN(NA) 
C 

IF  (IZ  .EQ.  0)  THEN 
CALL  DLSACG ( IXCRNT , CXPQNO , KXCRT , CIWO , 1 , CKLNO ) 
CALL  ESCFAR (NA, CKLNO , CKLN, CETHN, CEPHN) 
C 

C   If  we  calculate  in  X  and  Y  components  as  theta  approaches  0  or  pi, 
C   then  theta  component  is  cosin  phi  in  X  direction  plus  sine  phi  in  Y 
C   direction,  and  phi  component  is  negative  sine  phi  in  X  direction  plus 
C   cosin  phi  in  Y  direction. 

IF  (RTHE  .EQ.  0 . DO )   THEN 
CETHN=COS(RPHI) *CETHN+SIN(RPHI ) * CEPHN 
CEPHN= -SIN (RPHI ) *CETHN+COS (RPHI ) *CEPHN 

ELSEIF  (RTHE  . EQ .  PI)  THEN 
CETHN=- COS (RPHI) 'CETHN-SIN (RPHI ) * CEPHN 
CEPHN=-SIN(RPHI) 'CETHN+COS (RPHI ) *CEPHN 
END  IF 
C 

CESTH=CESTH+CETHN 
CESPH=CESPH+CEPHN 

IF  (NA  .NE.  0)  THEN 

IF  (IE  .EQ.  1)  THEN 
DO  2400  1=2, KXCRT, 2 
CKLNO (I)=-CKLN0 (I) 
2400   CONTINUE 

END  IF 

IF  (IM  .EQ.  1)  THEN 
DO  2500  1=1, KXCRT, 2 
CKLNO(I)=-CKLN0 (I) 
2500   CONTINUE 

END  IF 
NA1=-NA 
DN=NA1 
CALL  ESCFAR (NA1 , CKLNO , CKLN, CETHN, CEPHN) 

IF  (RTHE  .EQ.  0 . DO )   THEN 
CETHN=COS (RPHI ) *CETHN+SIN (RPHI ) *CEPHN 
CEPHN=-SIN (RPHI ) *CETHN+COS (RPHI ) *CEPHN 

ELSEIF  (RTHE  .EQ.  PI)  THEN 
CETHN=-COS (RPHI ) *CETHN-SIN (RPHI ) * CEPHN 
CEPHN=-SIN (RPHI ) *CETHN+COS (RPHI ) 'CEPHN 
END  IF 
C 

CESTH=CESTH+CETHN 
CESPH=CESPH+CEPHN 
END  IF 
END  IF 
C 

IF  (IZ  .EQ.  1)  THEN 
DO  3100  I=1,KCRNT 
DO  3000  J=1,KCRNT 

CXRPQ ( I , J) =CXPQN ( I , J) +CRPQ1 ( I , J) 
3  000     CONTINUE 
3100   CONTINUE 
DN=NA 

CALL  DLSACG (ICRNT, CXRPQ, KCRNT,CIW, l.CKLN) 
C 

IF  ((IK  .EQ.  1)  .AND.   (IR  . EQ .  45))  THEN 
CALL  KLCRNT(DN,CKLN,CIW) 
END  IF 
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CALL  ESCFAR (NA, CKLNO , CKLN, CETHN, CEPHN) 

IF  (RTHE  .EQ.  0 . DO )   THEN 
CETHN=COS(RPHI) *CETHN+SIN (RPHI ) *CEPHN 
CEPHN=-SIN(RPHI) *CETHN+COS (RPHI ) * CEPHN 

ELSEIF  (RTHE  .EQ.  PI)  THEN 
CETHN= -COS (RPHI) *CETHN-SIN (RPHI ) *CEPHN 
CEPHN=-SIN(RPHI) *CETHN+COS (RPHI ) *CEPHN 
END  IF 
C 

CESTH=CESTH+CETHN 
CESPH=CESPH+CEPHN 

IF  (NA  .NE.  0)  THEN 

IF  (NYSM  .EQ.  0)  THEN 
DO  3300  I=1,KCRNT 

DO  3200  J=1,KCRNT 
CXRPQ ( I , J) =CXPQN ( I , J ) +CRPQ2 ( I ,  J ) 
3200     CONTINUE 
3  3  00   CONTINUE 

ELSE 
DO  3  500  I=1,KCRNT 

DO  3400  J=1,KCRNT 
CXRPQ ( I , J) =CXPQN ( I , J ) +CRPQ1 ( I , J) 
3400     CONTINUE 
3  50  0   CONTINUE 

END  IF 
CALL  DLSACG (ICRNT, CXRPQ, KCRNT, CIW, 1 , CKLN) 

IF  (IE  .EQ.  1)  THEN 
DO  3600  1=2, KCRNT, 4 
11=1+1 

CKLN(I)=-CKLN(I) 
CKLN(I1)=-CKLN(I1) 
3  600   CONTINUE 

END  IF 

IF  (IM  .EQ.  1)  THEN 
DO  3700  1=1, KCRNT, 4 
12=1+3 

CKLN(I)=-CKLN(I) 
CKLN(I2)=-CKLN(I2) 
37  00   CONTINUE 

END  IF 
NA1=-NA 
DN=NA1 
C 

IF  ((IK  .EQ.  1)  .AND.  (IS  . EQ .  199))  THEN 
CALL  KLCRNT ( DN , CKLN , CIW ) 
END  IF 
C 

CALL  ESCFAR (NA1 , CKLNO , CKLN, CETHN, CEPHN) 

IF  (RTHE  .EQ.  0 . DO )   THEN 
CETHN=COS(RPHI) *CETHN+SIN (RPHI) *CEPHN 
CEPHN=- SIN (RPHI) *CETHN+COS (RPHI ) *CEPHN 

ELSEIF  (RTHE  . EQ .  PI)  THEN 
CETHN=-COS(RPHI) 'CETHN-SIN (RPHI ) *CEPHN 
CEPHN= - SIN ( RPHI ) *  CETHN+COS ( RPHI ) *  CEPHN 
END  IF 
C 

CESTH=CESTH+CETHN 
CESPH=CESPH+CEPHN 
END  IF 
END  IF 
C 

IF  (IK  .EQ.  1)  THEN 

IF  (IS  .EQ.  199)  THEN 
OPEN(31,FILE= 'khzz41.dz ' , STATUS =' UNKNOWN ' ) 
OPEN(3  2,FILE='lhzz41.dz ' , STATUS =' UNKNOWN ' ) 

END  IF 
DO  4200  IF=-30,30 

DO  4100  JZ=-30,30 
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CC 


AKPHPR=REAL ( CKPHP ( IF , JZ ) ) 
AKPHPI=IMAG(CKPHP(IF, JZ) ) 
AKZPR=REAL(CKZP(IF,  JZ)  ) 
AKZPI  =  IMAG(CKZP(IF,  JZ)  ) 
AKPHNR=REAL ( CKPHN ( IF , JZ ) ) 
AKPHNI=IMAG(CKPHN(IF, JZ) ) 
AKZNR=REAL ( CKZN ( IF , JZ ) ) 
AKZNI=IMAG(CKZN(IF, JZ) ) 
ALPHPR=REAL ( CLPHP ( IF , JZ )  ) 
ALPHPI  =  IMAG(CLPHP(IF, JZ)  ) 
ALZPR=REAL(CLZP(IF, JZ) ) 
ALZPI=IMAG(CLZP(IF, JZ) ) 
ALPHNR=REAL ( CLPHN ( IF , JZ ) ) 
ALPHNI  =  IMAG(CLPHN(IF, JZ)  ) 
ALZNR=REAL ( CLZN ( IF , JZ ) ) 
ALZNI  =  IMAG(CLZN(IF, JZ)  ) 


WRITE (31, *) 

AKPHPR, ' 

' ,AKPHPI 

WRITE (31, *) 

AKZPR, ' 

' ,AKZPI 

WRITE (32, *) 

ALPHPR, ' 

'  ,ALPHPI 

WRITE(32, *) 

ALZPR, ' 

' ,ALZPI 

WRITE (31,*) 

AKPHNR, ' 

' , AKPHNI 

WRITE (31, *) 

AKZNR, ' 

■ , AKZNI 

WRITE (32, *) 

ALPHNR, ' 

' ,ALPHNI 

WRITE (32, *) 

ALZNR, ' 

' ,ALZNI 

4100 

CONTINUE 

4200 

CONTINUE 
CLOSE (31) 
CLOSE (32) 
END  IF 
RETURN 
END 

SUBROUTINE  INCIDNT (NA, CIWO , CIW) 


C  This  subroutine  sets  up  the  incident  wave  on  an  object  which  is  coated 
C  with  anisotropic  material,  or  a  perfect  conductor. 

INCLUDE  'REALTP.INC 
INCLUDE  ' CMPXTP . INC ' 
INCLUDE  'MAINDM.INC 


DIMENSION  CIWO (KXCRT) ,CIW(KCRNT) , DJNS (KNDIM1+1 ) , DJPS (KQDIM1+2 ) 

COMMON  /INPUT2/  IE, IM, THETAI , THESINI , THECOSI , RTHEI 

COMMON  /INPUT4/  IZ , IK, IS, NYSM 

COMMON  /INPUT3/  RTHE, RDELT, RPHI , RDELP, NTHTAO, NPHI , THESIN, THECOS, 

►  RHPI 

COMMON  /RTHFTA/  DL1COSI , DL2SINI , DLL 

COMMON  /GCONST/  DKH, DKA, HH , HA, HSQ, DASQ, HSQN, DASQN, HHSQ, HDASQ, 

►  RAH,RAHSQ,DHA,DAH 

COMMON  /CRNTDM/  NMAX, MXNG , IQMAX, IQMAX1 , IQMAX2 , IXCRNT, ICRNT, MXQEG, 
•■  MXQOG 

NP=NA+1 

NP1=NP+1 

NM=NA-1 


IF  (IZ  .EQ.  0)  GO  TO  4000 
C 

C   This  part  computes  the  incident  field  on  an  anisotropic  object 
C 
C   Initialize  the  column  matrix  of  sum  current  on  the  anisotropic 

DO  100  IJ=1,KCRNT 

CIW(IJ)=CZERO 
100     CONTINUE 
C 

C   If  the  incident  angle  is  90  degree  or  0  degree 
IF  (RTHEI  .EQ.  ZERO)  GO  TO  3000 

CALL  DBSJNS(DL2SINI,  KNDIM1 ,  DJNS) 
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DJN=DJNS(NP) 

IF  (NA.  EQ.  0)  THEN 
DDJN=-DJNS(NP1) 

ELSE 
DDJN=HALF* (DJNS (NA) -DJNS (NP1) ) 

ENDIF 
C  Use  IMSL  library  to  compute  Bessel ' s  function  series. 
CALL  DBSJNS(DL1C0SI,IQMAX2,DJPS) 
DO  600  IP=0,KQDIM 
NIP=IP+1 
DNIP=NIP 
NIP1=NIP+1 
NIP2=NIP+2 
N11=NA+IP-1 
N12=NA+IP 
IE1=M0D(N11,4) 
IE2=MOD(N12,4) 
C 

CF1=CI2**IE1 
CF2=CI2**IE2 
IPW1=4*IP+1 
IPW2=IPW1+1 
IPW3=IPW2+1 
IPW4=IPW3+1 
C 

DJP=DJPS(NIP) 

DJP1=DJPS(NIP1) 

DJP2=DJPS(NIP2) 

IF  (IE  .EQ.  1)  THEN 
CEEPHI=CF1* (DJP+DJP2) *DDJN/DNIP 
CEHPHI=-TW0*NA*CF2*DJP1*DJN/DLL 
CEHZ=-CF2*THESINI* (DJP+DJP2) *DJN/DNIP 

ELSE 
CEEPHI=CZERO 
CEHPHI=CZERO 
CEHZ=CZERO 

END  IF 
C 

IF  (IM  .EQ.  1)  THEN 
CMEPHI=TW0*NA*CF2*DJP1*DJN/DLL 
CMEZ=CF2*THESINI* (DJP+DJP2) *DJN/DNIP 
CMHPHI=CF1* (DJP+DJP2 ) *DDJN/DNIP 

ELSE 
CMEPHI=CZERO 
CMEZ=CZERO 
CMHPHI=CZERO 

END  IF 
C 

CIW(IPW1)=TW0* (CEEPHI+*CMEPHI) 
CIW ( IPW2 ) =TWO*CMEZ 
CIW ( IPW3 ) =TWO* (CEHPHI+CMHPHI ) 
CIW ( IPW4 ) =TWO*CEHZ 
600    CONTINUE 

RETURN 
C 

3000   CALL  DBSJNS(DKH, IQMAX2,DJPS) 
DO  3200  IP=0,KQDIM 
JP=IP+2 
DJP=DJPS(JP) 
IEll=MOD(IP,4) 
IE12=MOD(IP+l,4) 
CF1=CI2**IE11 
CF2=CI2**IE12 
IPW1=4*IP+1 
IPW3=IPWl+2 
C 

IF  (IE  .EQ.  1)  THEN 
CEEPHI=CF1*DJP/DKH 
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CEHPHI=-CF2*DJP/DKH 

ELSE 
CEEPHI=CZERO 
CEHPHI=CZERO 
END  IF 
C 

IF  (IM  .EQ.  1)  THEN 
CMEPHI=CF2*DJP/DKH 
CMHPHI=CF1*DJP/DKH 

ELSE 
CMEPHI=CZERO 
CMHPHI=CZERO 
END  IF 
C 

CIW(IPWl) =TWO* (CEEPHI+CMEPHI) 
CIW(IPW3)=TW0* (CEHPHI+CMHPHI) 
3200   CONTINUE 

RETURN 
C 

C   This  part  compute  incident  fields  on  a  perfect  conductor 
C   Initialize  the  column  matrix  of  sum  current  on  the  conductor 
4000   DO  4200  IX=1,KXCRT 

CIW0 (IX)=CZERO 
4200    CONTINUE 
C 
C   If  the  incident  angle  is  90  degree  or  0  degree 

IF  (RTHEI  .EQ.  ZERO)  GO  TO  6000 
C 

CALL  DBSJNS(DL2SINI,  MXNG+1,  DJNS) 
DJN=DJNS(NP) 

IF  (NA  .EQ.  0)  THEN 
DDJN=-DJNS (NP1) 

ELSE 
DDJN=HALF* (DJNS(NA) -DJNS(NPl) ) 
ENDIF 
C 

CALL  DBSJNS(DLlCOSI,IQMAX2,DJPS) 
DO  4400  IP=0,KQDIM 
NIP=IP+1 
DNIP=NIP 
NIP1=NIP+1 
NIP2=NIP+2 
N11=NA+IP-1 
N12=NA+IP 
IEl=MOD(Nll, 4) 
IE2=MOD(N12, 4) 
C 

CF1=CI2**IE1 
CF2=CI2**IE2 
IPW1=2*IP+1 
IPW2=IPW1+1 
C 

DJP=DJPS(NIP) 
DJP1=DJPS (NIP1) 
DJP2=DJPS(NIP2) 

IF  (IE  .EQ.  1)  THEN 
CEEPHI=CF1* (DJP+DJP2) *DDJN/DNIP 

ELSE 
CEEPHI=CZERO 
END  IF 
C 

IF  (IM  .EQ.  1)  THEN 
CMEPHI=TWO*NA*CF2*DJPl*DJN/DLL 
CMEZ=CF2*THESINI* (DJP+DJP2) *DJN/DNIP 

ELSE 
CMEPHI=CZERO 
CMEZ=CZERO 
END  IF 
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CIWO (IPW1) =TWO* (CEEPHI+CMEPHI) 

CIWO ( IPW2 ) =TWO*CMEZ 
4400    CONTINUE 

RETURN 
C 
6000   CALL  DBS JNS(DKH, IQMAX2 , DJPS) 

DO  6200  IP=0,KQDIM 

JP=IP+2 

DJP=DJPS(JP) 

IE11=M0D(IP,4) 

IE12=MOD(IP+l,4) 

CF1=CI2**IE11 

CF2=CI2**IE12 

IPW1=2*IP+1 
C 

IF  (IE  .EQ.  1)  THEN 

CEEPHI=CF1*DJP/DKH 
ELSE 

CEEPHI=CZERO 
END  IF 
C 

IF  (IM  .EQ.  1)  THEN 

CMEPHI=CF2*DJP/DKH 
ELSE 

CMEPHI=CZERO 
END  IF 
C 

CIWO (IPW1)=TW0* (CEEPHI+CMEPHI) 
62  00   CONTINUE 

RETURN 

END 
C 

SUBROUTINE  XPQO 

Q*  **********  *  *  *  *  *  *  *  *  *  *  *  *  *  *****  *  *****  *  *  *  *  *************************  *  *  *  *  *-* 

C   This  subroutine  computes  the  matrix  XN(P,Q)  for  N  =  0  following  a 
C   call  to  XPQINI .  This  matrix  is  kept  in  the  common  block: 
C      COMMON  /XPQTMP/  CXPQN 

INCLUDE  ' REALTP . INC ' 

INCLUDE  ' CMPXTP . INC ' 

INCLUDE  ' MAINDM . INC ' 
C 

DIMENS I ON  CXPQN ( KCRNT , KCRNT ) , CXRPQ ( KCRNT , KCRNT ) 

DIMENSION  CXPQN0 (KXCRT , KXCRT) 

DIMENSION  CGNE ( 0 : MAXPEG , 0 :MAXPEG, KNDIM1  +  1 )  , 
+  CGNO ( 0 : MAXPOG , 0 :MAXPOG, KNDIM1+1 ) 

DIMENSION  CDGNE ( 0 : MAXPEG , 0: MAXPEG, KNDIM1  +  1)  , 
+  CDGNO ( 0 : MAXPOG , 0 : MAXPOG , KNDIM1+1 ) 

COMMON  /GPQTMP/  CGNE , CDGNE , CGNO , CDGNO 

COMMON  /XPQTMP/  CXPQN, CXRPQ, CXPQN0 

COMMON  /CRNTDM/  NMAX, MXNG , IQMAX, IQMAX1 , IQMAX2 , IXCRNT, ICRNT, MXQEG , 
+  MXQOG 

COMMON  /GCONST/  DKH , DKA, HH, HA, HSQ, DASQ, HSQN, DASQN, HHSQ, HDASQ, 
+  RAH,RAHSQ,DHA,DAH 

COMMON  /INPUT4/  IZ , IK, IS , NYSM 

SAVE  /XPQTMP/ , /GPQTMP/, /GCONST/ 
C 

IF  (N  .NE.  0)  THEN 

WRITE (*,*)  'Input  N  is  not  equal  to  0  in  XPQO.' 

WRITE!*,*)  'Execution  is  stopped.' 

STOP 

END  IF 
C 

IF  (IZ  .EQ.  0)  GO  TO  4000 
C 

C   Initialize  the  XN(P,Q)  matrix. 
C   Null  terms  will  be  skipped  later. 
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100 
200 
C 
C 


DO  200  IQ=1,KCRNT 

DO  100  IP=1,KCRNT 

CXPQNdP,  IQ)=CZERO 
CONTINUE 

CONTINUE 


Form  the  XN(P,Q) 
CF31=DKA*CI1 


matrix  for 


0  and  using  on  a  perfect  conductor 


DO  1300  IQE=0,MXQEG-1 

IQE1=IQE+1 

IQ=2*IQE 

IQX1=4*IQ+1 

IQX2=IQX1+1 

IQX3=IQX2+1 

IQX4=IQX3+1 

DQ1=IQ+1 

FQ22=DQ1/HHSQ 

DO  1100  IPE=0,MXQEG-1 

IPE1=IPE+1 

IP=2*IPE 

IPX1=4*IP+1 

IPX2=IPX1+1 

IPX3=IPX2+1 

IPX4=IPX3+1 

DP1=IP+1 

F41=DKH/DP1 

CF11=F41*CF31 

F51=QUAR*F41*DKA 

CXPQN (IPX1, IQX1) = (CGNEdPEl, IQE, 2) -CGNE (IPE, IQE, 2 ) ) *CF11 

CXPQN(IPX4, IQX1) = (CGNE(IPE1,IQE, 2) -CGNE (IPE, IQE, 2) +HA* ( 
+  CDGNE(IPE1,IQE,2) -CDGNE (IPE, IQE, 2 ) ) ) *F41 

CXPQN (IPX2,IQX2)=(TWO*FQ22*DPl*CGNO( IPE, IQE, D+HALF* ( 
+  CGNE (IPE, IQE1, 1)+CGNE(IPE1, IQE, 1) - 

+  CGNE(IPE1,IQE1, 1) -CGNE (IPE, IQE,1) ) ) *CF11 

CXPQN(IPX3 , IQX2)= (CDGNE (IPE, IQE, 1 ) +CDGNE (IPE1 , IQE1 , 1) - 
+  CDGNE (IPE1, IQE, 1) -CDGNE ( IPE, IQE1 , 1 ) ) *F51 

CXPQN ( IPX2 , IQX3 ) =-CXPQN ( IPX 4 , IQX1 ) 

CXPQN ( IPX3 , IQX3 ) =CXPQN ( IPX1 , IQX1 ) 

CXPQN (IPX1,IQX4)=-CXPQN(IPX3,IQX2) 

CXPQN (IPX4 , IQX4 ) =CXPQN ( IPX2 , IQX2 ) 
1100      CONTINUE 
1300   CONTINUE 

DO  2300  IQO=0,MXQOG-1 

IQ01=IQO+l 

IQ=2*IQO+l 

IQX1=4*IQ+1 

IQX2=IQX1+1 

IQX3=IQX2+1 

IQX4=IQX3+1 

DQ1=IQ+1 

FQ22=DQ1/HHSQ 

DO  2200  IPO=0,MXQOG-1 

IP01=IPO+l 

IP=2*IPO+l 

IPX1=4*IP+1 

IPX2=IPX1+1 

IPX3=IPX2+1 

IPX4=IPX3+1 

DP1=IP+1 

F41=DKH/DP1 

CF11=F41*CF31 

F51=QUAR*F41*DKA 

CXPQN(IPX1,IQX1)  =  (CGNO(IP01,IQO,2)-CGNO(IPO,IQO,2)  )  *CF11 

CXPQN ( IPX4 , IQX1 ) = (CGNO ( IPOl , IQO , 2 ) -CGNO ( IPO , IQO , 2 ) +HA* ( 
+  CDGNO(IP01,IQO,2)-CDGNO(IPO,IQO,2) ) ) *F41 

CXPQN (IPX2, IQX2)= (TW0*FQ22*DP1*CGNE ( IPOl , IQOl , 1) +HALF* ( 
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+  CGN0(IP0,IQ01,1)+CGN0(IP01,IQ0,1) - 

+  CGNO (IPOl, IQOl, 1) -CGN0(IP0,IQ0,1) ) ) *CF11 

CXPQN ( IPX3 , IQX2 ) = (CDGNO (IPO, 100, 1 ) +CDGNO ( IP01 , IQ01 , 1 ) - 
+  CDGN0(IP01,IQ0,1) -CDGN0(IP0,IQ01,1) ) *F51 

CXPQN ( IPX2 , IQX3 ) =-CXPQN ( IPX4 , IQX1 ) 

CXPQN (IPX3,IQX3) =CXPQN(IPX1, IQX1 ) 

CXPQN ( IPX1 , IQX4 ) =-CXPQN { IPX3 , IQX2 ) 

CXPQN ( IPX4 , IQX4 ) =CXPQN ( IPX2 , IQX2 ) 
2200      CONTINUE 
23  00   CONTINUE 

RETURN 
C   Initialize  the  XN(P,Q)  matrix. 
C  Null  terms  will  be  skipped  later. 
C 
4000   DO  4200  IQ=1,KXCRT 

DO  4100  IP=1,KXCRT 

CXPQN0 (IP, IQ) =CZERO 
4100      CONTINUE 
4200    CONTINUE 
C 
C   Form  the  XN(P,Q)  matrix  for  n  =  0  and  using  on  an  anisotropic  coat. 

CF31=DKA*CI1 
C 

DO  4500  IQE=0,MXQEG-1 

IQE1=IQE+1 

IQ=2*IQE 

IQX1=2*IQ+1 

IQX2=IQX1+1 

DQ1=IQ+1 

FQ22=DQ1/HHSQ 

DO  4400  IPE=0,MXQEG-1 

IPE1=IPE+1 

IP=2*IPE 

IPX1=2*IP+1 

IPX2=IPX1+1 

DP1=IP+1 

F41=DKH/DP1 

CF11=F41*CF31 

F51=QUAR*F41*DKA 

CXPQN0(IPX1,IQX1)=(CGNE(IPE1,IQE,2) -CGNE (IPE, IQE, 2 ) ) *CF11 

CXPQN0 (IPX2, IQX2) = (TWO*FQ22*DPl*CGNO ( IPE , IQE, 1 ) +HALF* ( 
+  CGNE (IPE, IQE1, 1)+CGNE(IPE1, IQE, 1) - 

+  CGNE(IPE1,IQE1,1) -CGNE (IPE, IQE, 1) ) ) *CF11 

4400      CONTINUE 
4  500   CONTINUE 

DO  5300  IQO=0,MXQOG-1 

IQ01=IQO*l 

IQ=2*IQO+l 

IQX1=2*IQ+1 

IQX2=IQX1+1 

DQ1=IQ+1 

FQ22=DQ1/HHSQ 

DO  5200  IPO=0,MXQOG-1 

IP01=IPO+l 

IP=2*IPO+l 

IPX1=2*IP+1 

IPX2=IPX1+1 

DP1=IP+1 

F41=DKH/DP1 

CF11=F41*CF31 

F51=QUAR*F41*DKA 

CXPQN0 (IPXl,IQXl)=(CGNO(IP01,IQO,2)-CGNO(IPO,IQO,2) ) *CF11 

CXPQN0 (IPX2, IQX2)= (TW0*FQ22*DP1*CGNE ( IPOl , IQOl , D+HALF* ( 
+  CGNO ( IPO , IQOl , 1 ) +CGNO (IPOl , IQO , 1 ) - 

+  CGNO (IPOl, IQOl, 1) -CGNO (IPO, IQO, 1) ) ) *CF11 

5200      CONTINUE 
53  00   CONTINUE 

RETURN 
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END 
C 

SUBROUTINE  XPQN(NIN) 

C   This  subroutine  calls  the  subroutine  GDN  to  update  G(P,Q,N)  for 

C   N  =  NIN+1,  then  forms  the  matrix  XN(P,Q)  for  N  =  NIN  >  0.  This  matrix 

C   is  kept  in  the  common  block: 

C      COMMON  /XPQTMP/  CXPQN 

p*********************************************************************** 

C=~==— =============================================================== 

C  WARNING: 

C  IT  IS  ASSUMED  THAT  XPQINI  AND  XPQN  FOR  N  FROM  1  TO  NIN-1  HAVE  BEEN 

C  CALLED  SO  THAT  G(P,Q,N)  AND  ITS  DERIVATIVE  FOR  N=NIN-1  AND  N=NIN 

C  HAVE  BEEN  STORED  IN  THE  COMMON  BLOCK  /GPQTMP/  WITH  PROPER  N-INDICES. 

0======================================================================= 

INCLUDE  ' REALTP . INC ' 

INCLUDE  ' CMPXTP . INC ' 

INCLUDE  ' MAINDM . INC ' 
C 

DIMENSION  CXPQN ( KCRNT , KCRNT ) , CXRPQ ( KCRNT , KCRNT ) 

DIMENSION  CXPQNO (KXCRT, KXCRT) 

DIMENSION  CGNE(0:MAXPEG,0:MAXPEG,KNDIM1+1) , 
+  CGNO ( 0 : MAXPOG , 0 :MAXPOG, KNDIM1+1) 

DIMENSION  CDGNE ( 0 : MAXPEG , 0 : MAXPEG , KNDIM1+1 ) , 
+  CDGNO(0:MAXPOG,0:MAXPOG,KNDIM1+1) 

COMMON  /GPQTMP/  CGNE , CDGNE , CGNO , CDGNO 

COMMON  /XPQTMP/  CXPQN, CXRPQ, CXPQNO 

COMMON  /CRNTDM/  NMAX, MXNG, IQMAX, IQMAX1 , IQMAX2 , IXCRNT, ICRNT,  MXQEG , 
+  MXQOG 

COMMON  /GCONST/  DKH, DKA, HH, HA, HSQ, DASQ, HSQN, DASQN, HHSQ, HDASQ, 
+  RAH , RAHSQ , DHA , DAH 

COMMON  /NCONST/  DN,DNH,N 

COMMON  /INPUT4/  IZ , IK, IS , NYSM 

SAVE  /XPQTMP/ , /GPQTMP/ , /GCONST/ 
C 

N=NIN 

NA=ABS (N) 

IF  (NA  .LT.  1)  THEN 

WRITE (*,*)  'Input  ABS(N)  is  less  than  1  in  XPQN.' 

WRITE (*,*)  'Execution  is  stopped.' 

STOP 

END  IF 
C 

N0I=NA+1 

NPI=N0I+1 

NMI=NA 
C 

DN0=NA 

DNP=N0I 

DNM=NA-1 
C 

IF  (IZ  .EQ.  0)  GO  TO  4000 
C 

C   Initialize  the  XN(P,Q)  matrix. 
C  Null  terms  will  be  skipped  later. 
C 

DO  800  IQ=1, KCRNT 

DO  700  IP=1, KCRNT 

CXPQN ( IP , IQ ) =CZERO 
7  00       CONTINUE 
800    CONTINUE 
C 


Form  the  XN(P,Q)  matrix. 
CF31=DKA*CI1 
F21=TWO*DN0 
FN11=F21*DN0/DASQ 

DO  13  00  IQE=0,MXQEG-1 


104 


IQE1=IQE+1 

IQ=2*IQE 

IQX1=4*IQ+1 

IQX2=IQX1+1 

IQX3=IQX2+1 

IQX4=IQX3+1 

DQ1=IQ+1 

FQ12=DN0*DQ1 

FQ22=DQ1/HHSQ 

DO    1100    IPE=0,MXQEG-1 

IPE1=IPE+1 

IP=2*IPE 

IPX1=4*IP+1 

IPX2=IPX1+1 

IPX3=IPX2+1 

IPX4=IPX3+1 

DP1=IP+1 

F41=HH/DP1 

CF11=F41*CF31 

F51=HALF*DKA*F41 

CXPQN(IPX1,IQX1)  =  (  (CGNEdPE,  IQE,N0I)  -CGNEdPEl,  IQE, NOI)  )  *FN11+ 
+  CGNEdPEl,  IQE,NMI)  -CGNE(IPE,  IQE, NMI )+ 

+  CGNE(IPE1,IQE,NPI)-CGNE(IPE,IQE,NPI) ) *CF11 

CXPQN(IPX4,IQX1)  =  (  (CGNEdPE,  IQE,  NMI )  -CGNE  ( IPE1 ,  IQE, NMI)  )  *DNM+ 
+  (CGNEdPEl,  IQE, NPI)  -CGNEdPE,  IQE, NPI)  )  *DNP+ 

+  HA*  (CDGNEdPEl,  IQE, NMI)  -CDGNEdPE,  IQE, NMI) + 

+  CDGNE  (IPE1,  IQE,  NPI) -CDGNEdPE,  IQE,  NPI)  )  )  *F41 

CXPQN(IPX2, IQX2) = (FQ22*DP1*CGN0 ( IPE, IQE,N0I) + 
+  CGNEdPE,  IQE1,N0I)+CGNE(IPE1,  IQE.NOI)  - 

+  CGNEdPEl,  IQE1.N0I)  -CGNEdPE,  IQE.NOI)  )  *CF11 

CXPQN ( IPX3 , IQX2 ) = ( CDGNE ( IPE , IQE , NO I ) +CDGNE ( IPE1 , IQE1 , NO I ) - 
+  CDGNEdPEl,  IQE,  NOI)  -CDGNE  (IPE,  IQE1 ,  NOI )  )  *F51 

CXPQN ( IPX2 , IQX3 ) =-CXPQN ( IPX4 , IQX1 ) 

CXPQN ( IPX3 , IQX3 ) =CXPQN ( IPX1 , IQX1 ) 

CXPQN (IPX1, IQX4)=-CXPQN(IPX3 , IQX2) 

CXPQN ( IPX4 , IQX4 ) =CXPQN ( IPX2 , IQX2 ) 
1100      CONTINUE 

DO  1200  IPO=0,MXQOG-1 

IP01=IP0+1 

IP=2*IP0+1 

IPX1=4*IP+1 

IPX2=IPX1+1 

IPX3=IPX2+1 

IPX4=IPX3+1 

DP1=IP+1 

F12=FQ12/DP1 

CXPQN(IPX2,IQX1)=F21*CGNE(IP01, IQE.NOI) 

CXPQN  (IPX3,IQX1)  =  (CGNEdPOl,  IQE,  NMI)  -CGNEdPOl,  IQE,  NPI)  )  *CF31 

CXPQN ( IPX1 , IQX2 )  =  (CGNO ( IPOl , IQE , NOI ) -CGNO (IPO, IQE , NOI )  )  *F12 

CXPQN (IPX1 , IQX3 ) =-CXPQN ( IPX3 , IQX1 ) 

CXPQN ( IPX4 , IQX3 ) =CXPQN (IPX2 , IQX1 ) 

CXPQN ( IPX3 , IQX4 ) =CXPQN ( IPX1 , IQX2 ) 
1200      CONTINUE 
13  00   CONTINUE 

DO  2300  IQO=0,MXQOG-1 

IQ01=IQO+l 

IQ=2*IQO+l 

IQX1=4*IQ+1 

IQX2=IQX1+1 

IQX3=IQX2+1 

IQX4=IQX3+1 

DQ1=IQ+1 

FQ12=DN0*DQ1 

FQ22=DQ1/HHSQ 

DO  2100  IPE=0,MXQEG-1 

IPE1=IPE+1 

IP=2*IPE 

IPX1=4*IP+1 
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IPX2=IPX1+1 
IPX3=IPX2+1 
IPX4=IPX3+1 
DP1=IP+1 
F12=FQ12/DP1 

CXPQN ( IPX2 , IQX1 ) =F2 1 *CGNO ( IPE , IQO , NO I ) 

CXPQN ( IPX3 , IQX1 )  =  ( CGNO ( I PE , IQO , NMI ) -CGNO ( I PE , IQO , NPI ) ) *  CF3 1 
CXPQN(IPX1,IQX2)  =  (CGNEdPEl,  IQOl, NO I)  -CGNEdPE,  IQOl, NO I)  )  *F12 
CXPQN ( IPX1 , IQX3 ) =-CXPQN ( IPX3 , IQX1 ) 
CXPQN ( IPX4 , IQX3 ) =CXPQN ( IPX2 , IQX1 ) 
CXPQN ( IPX3 , IQX4 ) =CXPQN ( IPX1 , IQX2  ) 
2100      CONTINUE 

DO  2200  IPO=0,MXQOG-1 
IP01=IPO+l 
IP=2*IPO+l 
IPX1=4*IP+1 
IPX2=IPX1+1 
IPX3=IPX2+1 
IPX4=IPX3+1 
DP1=IP+1 
F41=HH/DP1 
CF11=F41*CF31 
F51=HALF*DKA*F41 

CXPQN(IPX1,IQX1)=( ( CGNO ( IPO, IQO, NO I) -CGNO (IPOl, IQO, NO I) ) *FN11+ 
+  CGNO (IPOl, IQO, NMI) -CGNO (IPO, IQO, NMI )+ 

+  CGNOdPOl,  IQO, NPI)  -CGNOdPO,  IQO, NPI)  )  *CF11 

CXPQN  (IPX4,  IQX1)  =  (  (CGNOdPO,  IQO,  NMI)  -CGNO  (IPOl,  IQO,  NMI)  )  *DNM+ 
+  (CGNOdPOl,  IQO, NPI)  -CGNOdPO,  IQO, NPI)  )  *DNP+ 

+  HA*  (CDGNOdPOl,  IQO,  NMI)  -CDGNO  ( IPO,  IQO,  NMI )  + 

+  CDGNO (IPOl, IQO, NPI) -CDGNO (IPO, IQO, NPI) ) ) *F41 

CXPQN (IPX2, IQX2) = (FQ22*DP1*CGNE (IPOl , IQOl ,N0I ) + 
+  CGNO ( IPO, IQOl , NOI ) +CGNO (IPOl , IQO, NOI ) - 

+  CGNOdPOl,  IQOl, NOI)  -CGNOdPO,  IQO, NOI)  )  *CF11 

CXPQN(IPX3,IQX2)  =  (CDGNO (IPO, IQO, NOI ) +CDGNO (IPOl , IQOl,  NO I) - 
+  CDGNO (IPOl, IQO, NO I) -CDGNO (IPO, IQOl, NO I) ) *F51 

CXPQN ( IPX2 , IQX3 ) =-CXPQN ( IPX4 , IQX1 ) 
CXPQN (IPX3 , IQX3 ) =CXPQN ( IPX1 , IQX1 ) 
CXPQN ( IPX1 , IQX4 ) =-CXPQN ( IPX3 , IQX2 ) 
CXPQN ( IPX4 , IQX4 ) =CXPQN ( I PX2 , IQX2 ) 
220  0      CONTINUE 
2300   CONTINUE 

RETURN 
C   Initialize  the  XN(P,Q)  matrix. 
C  Null  terms  will  be  skipped  later. 
C 
4000   DO  4800  IQ=1,KXCRT 

DO  4700  IP=1,KXCRT 
CXPQN0 (IP, IQ) =CZERO 
47  0  0      CONTINUE 
480  0   CONTINUE 
C 

C   Form  the  XN(P,Q)  matrix. 
CF31=DKA*CI1 
F21=TWO*DN0 
FN11=F21*DN0/DASQ 
C 

DO  5300  IQE=0,MXQEG-1 

IQE1=IQE+1 

IQ=2*IQE 

IQX1=2*IQ+1 

IQX2=IQX1+1 

DQ1=IQ+1 

FQ12=DN0*DQ1 

FQ22=DQ1/HHSQ 

DO  5100  IPE=0,MXQEG-1 
IPE1=IPE+1 
IP=2*IPE 
IPX1=2*IP+1 
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5100 


5200 
5300 


6100 


6200 
6300 


IPX2=IPX1+1 

DP1=IP+1 

F41=HH/DP1 

CF11=F41*CF31 

F51=HALF*DKA*F41 

CXPQNO (IPX1,IQX1)=( (CGNE(IPE,IQE,N0I)-CGNE(IPE1,IQE,N0I) ) *FN11+ 
k  CGNE ( I PE1 , IQE , NMI ) -CGNE ( I PE , IQE , NMI ) + 

k  CGNE  (IPE1,  IQE,  NPI)  -CGNE (IPE,  IQE,NPI)  )  *CF11 

CXPQNO (IPX2, IQX2) = (FQ22*DP1*CGN0 (IPE, IQE , NOI ) + 
k  CGNE (IPE , IQE1 , NOI ) +CGNE ( IPE1 , IQE, NOI ) - 

k  CGNEdPEl,  IQE1, NOI)  -CGNE(IPE,  IQE, NOI)  )  *CF11 

CONTINUE 
DO  5200  IPO=0,MXQOG-1 

IP01=IPO+l 

IP=2*IPO+l 

IPX1=2*IP+1 

IPX2=IPX1+1 

DP1=IP+1 

F12=FQ12/DP1 

CXPQNO ( IPX2 , IQX1 ) =F21*CGNE ( IP01 , IQE , NOI ) 

CXPQNO ( IPX1 , IQX2 ) = ( CGNO ( I POl , IQE , NO I ) -CGNO ( IPO , IQE , NOI ) ) *F12 
CONTINUE 

CONTINUE 

DO  6300  IQO=0,MXQOG-1 

IQ01=IQO+l 

IQ=2*IQO+l 

IQX1=2*IQ+1 

IQX2=IQX1+1 

DQ1=IQ+1 

FQ12=DN0*DQ1 

FQ22=DQ1/HHSQ 

DO  6100  IPE=0,MXQEG-1 

IPE1=IPE+1 

IP=2*IPE 

IPX1=2*IP+1 

IPX2=IPX1+1 

DP1=IP+1 

F12=FQ12/DP1 

CXPQNO (IPX2,IQXl)=F21*CGNO(IPE,IQO,N0I) 

CXPQNO ( IPX1 , IQX2 ) = (CGNE ( IPE1 , IQOl ,N0I ) -CGNE (IPE, IQOl , NOI ) ) *F12 
CONTINUE 
DO  6200  IPO=0,MXQOG-1 

IP01=IPO+l 

IP=2*IPO+l 

IPX1=2*IP+1 

IPX2=IPX1+1 

DP1=IP+1 

F41=HH/DP1 

CF11=F41*CF31 

F51=HALF*DKA*F41 

CXPQNO (IPX1, IQX1)=( (CGNO ( IPO, IQO, NOI) -CGNO (IP01,IQO, NOI) ) *FN11+ 
h  CGNOdPOl,  IQO,  NMI)  -CGNO  ( IPO,  IQO,  NMI )  + 

^  CGNOdPOl,  IQO,  NPI)  -CGNO  (IPO,  IQO,  NPI)  )  *CF11 

CXPQNO (IPX2,IQX2)=(FQ22*DP1*CGNE(IP01, IQOl, NOI )+ 
^  CGNO (IPO, IQOl,N0I)+CGNO(IPOl, IQO,  NOI)  - 

*  CGNOdPOl,  IQOl, NOI)  -CGNO(IPO,  IQO, NOI)  )  *CF11 

CONTINUE 

CONTINUE 

RETURN 

END 


SUBROUTINE  ESCFAR (NIN, CKLN0 , CKLN, CETHN, CEPHN) 
This  subroutine  computes  the  scattered  fields  in  far  zone 


r*********i 


INCLUDE 
INCLUDE 
INCLUDE 


' REALTP . INC ' 
' CMPXTP . INC ' 
'MAINDM.INC 
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DIMENSION  CKLNO(KXCRT) 

DIMENSION  CKLN(KCRNT) , DJNS (KNDIM1 ) , DJPS (KQDIM1+2 ) 

COMMON  /CRNTDM/  NMAX, MXNG, IQMAX, IQMAX1 , IQMAX2 , IXCRNT, ICRNT,  MXQEG  , 
h  MXQOG 

COMMON  /GCONST/  DKH, DKA, HH, HA, HSQ, DASQ, HSQN, DASQN, HHSQ,  HDASQ, 
h  RAH,RAHSQ,DHA,DAH 

COMMON  /INPUT3/  RTHE, RDELT, RPHI , RDELP , NTHTAO , NPHI , THESIN, THECOS, 
h  RHPI 

COMMON  /INPUT4/  IZ, IK, IS,NYSM 

N=NIN 

DN1=N 

NA=ABS (NIN) 

NP=NA+1 

NP1=NP+1 

PHI1=RPHI*N 

CI2N=CI2**N 

CENP=CONE*COS(PHIl)+CIl*SIN(PHIl) 

CETHN=CZERO 

CEPHN=CZERO 

IF  ((RTHE  .EQ.  0.)  .OR.  (RTHE  . EQ .  PI))  GO  TO  2000 

TCOT=THECOS /THESIN 

DL2SIN=DKA*THESIN 

DLlCOS=DKH*THECOS 

DNL2=DNl*TCOT/DKA 

CALL  DBSJNS(DL2SIN,  KNDIM1 ,  DJNS) 

DJN=DJNS(NP) 

IF  (N  .EQ.  0)  THEN 
DDJN=-DJNS(NP1) 

ELSE 
DDJN=HALF* (DJNS (NA) -DJNS (NP1) ) 

ENDIF 

KN=NA/2 
KNP=NP/2 

IF  ( (N  .LT.  0)   .AND.   (KNP  . GT .  KN) )  THEN 
DJN=-DJN 
DDJN=-DDJN 
ENDIF 

CALL  DBSJNSfDLlCOS,  KQDIM1+2,  DJPS) 


IF  (IZ  .EQ.  0)  THEN 
DO  100  IP=0, IQMAX 
INP=IP+1 
NP2=INP+2 
N11=IP 
N12=IP+2 
IEl=MOD(Nll, 4) 
CF1=CI2**IE1 
IPW1=2*IP+1 
IPW2=IPW1+1 
DJP=DJPS(INP) 
DJP2=DJPS(NP2) 
DJPH=HALF* (DJP+DJP2 ) 

CETHN=CETHN+CI2N*DHA*CENP*CI1*CF1*DJN* (DNL2*CKLN0 (IPW1) *DJP 
+       -THESIN*CKLN0 (IPW2)*DJPH) 
CEPHN=CEPHN-CI2N*DHA*CENP*DDJN*CKLN0 (IPW1) *CF1*DJP 
100    CONTINUE 
ELSE 
DO  200  IP=0, IQMAX 
INP=IP+1 
NP2=INP+2 
N11=IP 
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N12=IP+2 
IE1=M0D(N11,4) 
CF1=CI2**IE1 
IPW1=4*IP+1 
IPW2=IPW1+1 
IPW3=IPW2+1 
IPW4=IPW3+1 
DJP=DJPS(INP) 
DJP2=DJPS(NP2) 
DJPH=HALF* (DJP+DJP2  ) 

CETHN=CETHN+CI2N*DHA*CENP* ( (CI1*DJN*DNL2*CKLN ( IPW1 ) -DDJN 
+       *CKLN(IPW3) ) *DJP-CI1*DJN*THESIN*CKLN(IPW2) *DJPH) *CF1 

CEPHN=CEPHN-CI2N*DHA*CENP* ( (DDJN*CKLN ( IPW1 ) +CI1*DJN*DNL2 
+       *CKLiN(IPW3) ) *DJP-CI1*DJN*THESIN*CKLN(IPW4) *DJPH) *CF1 
2  00    CONTINUE 
END  IF 
RETURN 
C 

2000      IF  (NA  .NE.  1)  THEN 
GO  TO  3000 

END  IF 
CALL  DBSJNS(DKH,KQDIM1,DJPS) 
C 

IF  (IZ  .EQ.  0)  THEN 
DO  2100  IP=0,IQMAX 
INP=IP+1 
IPW1=2*IP+1 
JP=MOD(IP,4) 
CF1=CI1**JP 
CF2=CI2**JP 
C 

DJP=DJPS(INP) 

IF  (RTHE  .EQ.  0.)  THEN 
IF  (N  .LT.  0)  THEN 
CETHN=CETHN-DAH*DJP*CF2*CKLN0 (IPW1) 
CEPHN=CEPHN+DAH*DJP*CF2*CI1*CKLN0 (IPW1) 

ELSE 
CETHN=CETHN+DAH*DJP*CF2*CKLN0 (IPW1) 
CEPHN=CEPHN+DAH*DJP*CF2*CI1*CKLN0(IPW1) 
END  IF 
END  IF 

IF  (RTHE  .EQ.  PI)  THEN 
IF  (N  .LT.  0)  THEN 
CETHN=CETHN-DAH*DJP*CF1*CKLN0 ( IPW1 ) 
CEPHN=CEPHN+DAH*DJP*CF1*CI1*CKLN0 (IPW1) 

ELSE 
CETHN=CETHN+DAH*DJP*CF1*CKLN0 (IPW1) 
CEPHN=CEPHN+DAH*DJP*CF1*CI1*CKLN0 (IPW1) 
END  IF 
END  IF 
2100   CONTINUE 
ELSE 
DO  2200  IP=0,IQMAX 
INP=IP+1 
IPW1=4*IP+1 
IPW3=IPWl+2 
JP=MOD(IP,4) 
CF1=CI1**JP 
CF2=CI2**JP 
C 

DJP=DJPS(INP) 

IF  (RTHE  .EQ.  0.)  THEN 
IF  (N  .LT.  0)  THEN 
CETHN=CETHN-DAH»DJP*CF2* (CKLN(IPWl) -CI1*CKLN (IPW3  )  ) 
CEPHN=CEPHN+DAH*DJP*CF2  * (CI1*CKLN ( IPW1 ) +CKLN ( IPW3 )  ) 

ELSE 
CETHN=CETHN+DAH*DJP*CF2* (CKLN ( IPW1 ) +CI1*CKLN ( IPW3  )  ) 
CEPHN=CEPHN+DAH*DJP*CF2* (CI1*CKLN ( IPW1 ) -CKLN(IPW3) ) 
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END  IF 
END  IF 

IF  (RTHE  .EQ.  PI)  THEN 
IF  (N  .LT.  0)  THEN 
CETHN=CETHN-DAH*DJP*CF1* (CKLN (IPW1 ) +CI1*CKLN (IPW3  )  ) 
CEPHN=CEPHN+DAH*DJP*CF1* (CI1*CKLN (IPW1 ) -CKLN(IPW3) ) 

ELSE 
CETHN=CETHN+DAH*DJP*CF1* (CKLN (IPW1) -CI1*CKLN (IPW3 ) ) 
CEPHN=CEPHN+DAH*DJP*CF1* (CI1*CKLN (IPW1 ) +CKLN (IPW3 ) ) 
END  IF 
END  IF 
2200   CONTINUE 
END  IF 
3000   RETURN 
END 


SUBROUTINE  KLCRNT (DN, CKLN, CIW) 
p************************************************************************** 

C   This  subroutine  computes  equivalent  currents  K  and  L  on  inner  and  outer 
C  surfaces  respectively. 

INCLUDE  ' REALTP . INC ' 

INCLUDE  'CMPXTP.INC 

INCLUDE  ' MAINDM . INC ' 
C 

DIMENSION  CKLN ( KCRNT ), CIW (KCRNT) 

DIMENSION  CKLNP ( KCRNT ) , CKLNN ( KCRNT) 

DIMENSION  CRPQ31 (KCRNT, KCRNT) , CRPQ3 2 (KCRNT, KCRNT) 

DIMENSION  CKPHP(61,61) ,CKZP(61, 61) ,CLPHP(61, 61) , 
+  CLZP(61,61),  CKPHN(61, 61) ,CKZN(61, 61) , 

+  CLPHN(61, 61) ,CLZN(61, 61) 

COMMON  /CKLMTX/  CKPHP , CKZP , CLPHP, CLZP, CKPHN, CKZN, CLPHN, CLZN 

COMMON  /XPQTMP2/  CRPQ3 1 , CRPQ3  2 

COMMON  /CRNTDM/  NMAX, MXNG, IQMAX, IQMAX1 , IQMAX2 , IXCRNT , ICRNT  ,  MXQEG  , 
+  MXQOG 

COMMON  /INPUT3/  RTHE, RDELT, RPHI , RDELP , NTHTAO , NPHI , THESIN, THECOS, 
+  RHPI 

COMMON  /INPUT4/  IZ , IK, IS , NYSM 
C 

DO  400  IX=1, ICRNT 

CKLNP ( IX )=CZERO 

CKLNN (IX) =CZERO 
400     CONTINUE 

PHI1=DN*RPHI 

CPHI=CONE*COS(PHIl)+CIl*SIN(PHIl) 
C 

DO  600  IX=1, ICRNT 
DO  500  IY=1, ICRNT 

CKLNP(IX)=CKLNP(IX)+CRPQ31 (IX, IY) *CKLN(IY) 
500       CONTINUE 
600    CONTINUE 
C 

DO  800  IX=1, ICRNT 
DO  700  IY=1, ICRNT 

CKLNN(IX)=CKLNN(IX)+CRPQ32(IX,IY)*CKLN(IY) 
7  00       CONTINUE 
800    CONTINUE 
C 

DO  1600  IF=-30,30 

DF  =  IF 

DPHI=DF*PI/32.D0 

DNPHI =DN*DPHI 

CPHI=CONE*DCOS (DNPHI ) +CI1*DSIN (DNPHI ) 
DO    1500    JZ=-30,30 

DJZ=JZ 

DZ=DJZ/32 .DO 

DV=DACOS(DZ) 

SINV=DSIN(DV) 


110 


CKPHPdF,  JZ)  =CZERO 

CKZP(IF, JZ)=CZERO 

CLPHPdF,  JZ)=CZERO 

CLZPdF, JZ)=CZERO 

CKPHN ( IF , JZ ) =CZERO 

CKZN(IF, JZ)=CZERO 

CLPHN ( IF , JZ ) =CZERO 

CLZNdF,  JZ)=CZERO 
C 

DO  1400  IP=0,IQMAX 

P=IP 

P1=IP+1 

DPV=P*DV 

DP1V=P1*DV 

DCSPV=DCOS(DPV) 

SINP1V=DSIN(DP1V) 

IPX1=IP*4+1 

IPX2=IPX1+1 

IPX3=IPX2+1 

IPX4=IPX3+1 

CKPHP ( IF , JZ ) =CKPHP (IF , JZ ) +CPHI * ( CKLNP (IPX1 ) +CIW ( IPX4 ) ) 
+  *DCSPV/PI/SINV 

CKZP(IF, JZ)=CKZP(IF, JZ)+CPHI* (CKLNP (IPX2 ) -CIW(IPX3 ) ) *SINP1V/PI 

CLPHPdF,  JZ)=CLPHP  (IF,  JZ)+CPHI*  (CKLNP  (IPX3)  -CIW(IPX2)  ) 
+  *DCSPV/PI/SINV 

CLZPdF,  JZ)=CLZP(IF,  JZ)+CPHI*  (CKLNP  (IPX4  ) +CIW(IPX1)  )  *SINP1V/PI 

CKPHN (IF, JZ)=CKPHN(IF, JZ)+CPHI* (CKLNN(IPXl) -CIW(IPX4) ) 
+  *DCSPV/PI/SINV 

CKZNdF,  JZ)=CKZN(IF,  JZ)+CPHI*  (CKLNN(IPX2  )  +CIW(IPX3  )  )  *SINP1V/PI 

CLPHN ( IF , JZ ) =CLPHN ( IF , JZ ) +CPHI * ( CKLNN ( IPX3 ) +CIW ( IPX2 ) ) 
+  *DCSPV/PI/SINV 

CLZNdF,  JZ)=CLZN  (IF,  JZ)+CPHI*  (CKLNN  (IPX4)  -CIW(IPXl)  )  *SINP1V/PI 
1400         CONTINUE 
1500      CONTINUE 
1600   CONTINUE 

RETURN 

END 


J.         SUBROUTINE  RCSPAREA 


SUBROUTINE  RCSPAREA (CESTH, CESPH) 

C   This  subroutine  computes  cross  section  per  projected  area  in  all  direction. 
C 

INCLUDE  ' REALTP . INC ' 
INCLUDE  ' CMPXTP . INC ' 
C 

REAL* 4  ARCSPPA1 , ARCSPPA2 , APHASE1 , APHASE2 

COMMON  /GCONST/  DKH, DKA, HH , HA, HSQ, DASQ, HSQN, DASQN, HHSQ, HDASQ, 
+  RAH,RAHSQ,DHA,DAH 

COMMON  /INPUT3/  RTHE, RDELT , RPHI , RDELP, NTHTAO, NPHI , THESIN, THECOS, 
+  RHPI 

C 

ESTH=ABS (CESTH) 
ESPH=ABS (CESPH) 
ESTHSQ=ESTH*ESTH 
ESPHSQ=ESPH*ESPH 

IF  (RTHE  .EQ.  RHPI)  THEN 
RCSPPA1=PI*ESTHSQ/DKH/DKA 
RCSPPA2=PI*ESPHSQ/DKH/DKA 

ELSE  IF  ((RTHE  . EQ .  ZERO)  .OR.  (RTHE  . EQ .  PI))  THEN 
RCSPPA1=4.D0*ESTHSQ/DASQ 
RCSPPA2=4.D0*ESPHSQ/DASQ 
ELSE 


111 


AP=ABS(QUAR*DASQ*THECOS)+ABS(DKH*DKA*THESIN/PI) 

RCSPPA1=ESTHSQ/AP 

RCSPPA2=ESPHSQ/AP 

END  IF 
ARCSPPA1 =RCSPPA1 
ARCSPPA2  =RCSPPA2 

ER1=REAL(CESTH) 

EI1=IMAG(CESTH) 

ER2=REAL(CESPH) 

EI2=IMAG(CESPH) 

PHASE1=DATAN2 (EI1,ER1) 

PHASE2=DATAN2 (EI2.ER2) 

APHASE1=PHASE1 

APHASE2=PHASE2 

WRITE  (21,*)  ARCSPPA1,APHASE1,ARCSPPA2,APHASE2 

RETURN 
END 
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